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(a)  The  Generation  and  Evolution  of  Unsteady  Vorticity  in  a  Model  of  a  Solid 
Rocket  Engine  Chamber 

(b)  Nonlinear  Vorticity  Generation  by  Acoustic  Wave  Interaction  with  an  In¬ 
jected  Gas  Velocity  Field  in  a  Cylinder 

(c)  Unsteady  Vorticity  Generation  and  Evolution  in  a  Model  of  a  Solid  Rocket 
Engine  Chamber 


1.  Abstract 


A  new  formulation  for  chamber  flow  dynamics  in  a  model  of  solid  rocket  engine 
shows  that  vorticity  generation  and  convection  are  prominent  physical  features  of 
the  flow  field.  Analytical  and  fully  computational  methods  are  employed  to  describe 
a  basically  inviscid  interaction  between  acoustic  disturbances  arising  from  specified 
boundary  disturbances  and  a  sidewall  injected  flow  field  which  simulates  propellant 
burning.  The  mathematical  model,  based  on  the  Navier-Stokes  equations,  is  devel¬ 
oped  in  terms  of  an  initial  value  problem  in  order  to  describe  the  complete,  natural 
chamber  flow  evolution  arising  from  boundary  driven  disturbances.  The  approach 
is  analogous  to  a  direct  numerical  simulation,  although  cont  porary  perturbation 
methods  are  employed  to  extract  specific  spatial  and  temporal  scales  from  the  equa¬ 
tions  and  boundary  conditions.  The  results  show  that  large  unsteady  vorticity  is 
created  at  the  injected  surface  (sidewall)  and  convects  into  the  cylinder  with  the  ra¬ 
dial  component  of  the  injection  flow  velocity.  Eventually,  the  entire  chamber  is  filled 
with  an  intense,  rotational  flow  field  containing  relatively  large  velocity  gradients. 
The  presence  of  significant  distributed  vorticity  calls  into  question  the  predictions  of 
traditional  acoustic  stability  theory,  based  on  a  fundamentally  irrotational  concept. 
The  flow  field  predicted  by  the  present  theory  includes  an  essentially  acoustic  pres¬ 
sure  distribution,  and  a  velocity  distribution  composed  of  a  weakly  rotational  steady 
component,  an  irrotational  acoustic  component  and  a  strongly  rotational,  weakly  non¬ 
linear  and  viscous  component.  The  high  Reynolds  number,  low  Mach  number  flow 
field  has  many  similarities  to  a  conventionally  defined  internal  turbulent  flow. 
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2.  Project  Objectives,  Status 
and  Accomplishments 


Our  research  program  is  concerned  with  modelling  transient  flow  dynamics  in 
a  solid  rocket  engine  chamber.  Gaseous  mass  addition  from  propellant  burning  is 
simulated  by  fluid  injection  through  the  internal  wall  of  a  cylinder  with  one  closed  end. 
The  impact  of  spatially  distributed  burning  rate  transients  and/or  igniter  processes  on 
the  evolving  chamber  flow  is  assessed  from  the  perspective  of  an  initial  value  problem. 
In  particular,  chamber  flow  disturbances  are  created  by  specifically  defined,  spatially 
distributed,  time  dependent  variations  of  injection  velocity  or  pressure  on  boundary 
surfaces.  Solution  development  is  based  on  analytical  and  computational  methods, 
using  a  route  analogous  to  direct  numerical  simulation  (DNS),  where  the  entire  natural 
evolution  of  the  chamber  dynamics  is  described.  This  paradigm  shift  from  more 
familiar  acoustic  stability  theory  has  enabled  us  to  recognize  the  presence  of  new  and 
pertinent  physical  phenomena  in  the  chamber  flow.  In  particular,  we  demonstrate  that 
vorticity  is  generated  near  the  injection  (burning  propellant)  surface  and  transported 
into  the  entire  chamber  by  the  injected  flow  field.  The  presence  of  significant  levels 
of  time-dependent  vorticity  in  the  chamber  has  three  important  consequences.  First, 
traditional  acoustic  stability  theory,  involving  a  primarily  irrotational  disturbance 
flow  field,  cannot  describe  the  rotational  fields  found  in  our  analysis.  Secondly,  the 
shear  stress  near  the  propellant  surface  is  relatively  large  (compared  with  that  for  a 
time  averaged  mean  flow)  and  may  be  an  important  source  of  scouring  of  the  ^tcposed, 
degraded  propellant  surface.  Finally,  the  time-dependent  rotational  flow  field  in  our 
high  Reynolds  number,  low  Mach  number  chamber  has  similarities  to  turbulent  flow 
as  traditionally  conceived. 

During  the  project  period,  our  own  research  progress  enabled  us  to  create  very 
specific  programmatic  objectives: 

1.  Employ  the  complete  equations  of  motion  in  an  initial  value  analysis  to  study 
the  evolution  of  chamber  flow  transients  initiated  by  specified  boundary  driven 
disturbances. 

2.  Use  perturbation  methods,  based  on  viable  parameter  extremes  (large  Reynolds 
number,  small  Mach  number),  to  prove  that;  (1)  unsteady  vorticity  is  created 
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on  the  injection  surface  by  an  in  viscid  interaction  between  acoustic  waves  and 
the  injected  velocity  field,  and  (2)  the  vorticity  is  convected  into  the  chamber 
by  the  radial  component  of  the  injected  velocity. 

3.  Demonstrate  the  impact  of  relatively  weak  nonlinearity  and  viscous  effects  on 
the  spatial  distribution  of  the  transient  vorticity. 

4.  Relate  the  present  modelling  approach  for  solid  rocket  engine  chamber  dynamics 
to  the  traditional  acoustic  stability  theory. 

The  present  modelling  effort  is  used  to  demonstrate  how  and  where  unsteady 
vorticity  is  created  by  an  interaction  between  axially  propagating  acoustic  waves 
and  a  flow  field  caused  by  steady  sidewall  m£iss  injection  in  a  finite  length  cylinder. 
Additionally,  the  formulation  and  analysis  describes  the  transport  and  time-history  of 
the  spatial  distribution  of  vorticity  within  the  cylinder.  Finally,  the  results  show  that 
transient  rotational  flow  effects  are  crucial  to  the  evolution  and  stability  of  internal 
fluid  dynamics  when  the  characteristic  cylinder  Reynolds  number  (Re)  and  Mach 
number  (M)  are  very  large  and  small,  respectively. 

The  mathematical  model  first  describes  a  steady  flow  field  in  a  cylinder  with  one 
closed  end,  induced  by  strong  normal,  steady  sidewall  injection.  In  particular,  the 
injected  flow  is  characterized  by  a  radial  speed  much  larger  than  the  ratio  of  the  char¬ 
acteristic  axial  speed  to  Following  Taylor  (1956),  Culick  (1966)  shows  that 

the  steady  flow  is  described  in  a  first  approximation  by  an  inviscid,  rotational  equa¬ 
tion  system.  Axial  acoustic  disturbances  of  0[M)  are  created  by  either  a  prescribed 
time-dependent  axial  velocity  variation  at  the  closed  end  of  the  cylinder  or  a  similar 
pressure  variation  on  an  exit  boundary  upstream  of  the  chamber  nozzle.  They  prop¬ 
agate  through  the  btisic  inviscid  shear  flow  field,  and  perhaps  unexpectedly,  create 
significant  unsteady  vorticity  at  the  surface  of  the  porous  cylinder.  The  radial  com¬ 
ponent  of  the  steady  injected  flow  field  carries  the  vorticity  into  the  entire  cylinder. 
Earlier  related  work  (Wang  and  Kassoy,  1992  a,b)  demonstrates  that  refraction  of  the 
primary  axial  waves  can  generate  oblique  and  even  transverse  wave  disturbances  that 
propagate  through  the  shear  flow. 

Predictions  of  flow  dynamics  in  closely  related  systems  have  been  obtained  in 
the  past  from  traxlitional  acoustic  stability  theory  ^ls  formulated  and  described  by 
Grad  (1949).  Culick  and  co-workers,  (See  Culick  (1990)  for  a  comprehensive  review) 
in  particular,  have  developed  since  the  mid-1970’s  an  increasingly  sophisticated  lin¬ 
ear  and  weakly  nonlinear  stability  theory  to  describe  disturbance  behavior  observed 
in  laboratory  experiments  and  practical  systems  like  solid  rocket  engine  chambers. 
Comparisons  with  experiments  appear  to  be  reasonably  good,  although  the  theory  «is 
formulated  cannot  describe  the  generation  and  evolution  of  vorticity. 
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Recent  experimental  (Brown  et  al.  1986a, b,  Brown  and  ShaefFer,  1992)  and  com¬ 
putational  (Vuillot  and  Avalon, 1991;  Vuillot,1991)  studies  include  definitive  evidence 
of  significant  transient  rotational  flow  structures  distributed  throughout  the  injected 
cylinder  flow.  In  addition.  Price  and  Flandro  (1993)  and  Flandro  and  Roach  (1992) 
have  formulated  approximate  mathematical  models  for  describing  the  generation  and 
evolution  of  vorticity  in  an  injected  cylinder  flow  field.  These  contributions  provide 
the  motivation  for  developing  a  comprehensive  flow  dynamics  model  which  includes 
both  the  irrotational  acoustical  phenomena  of  traditional  stability  theory  and  impor¬ 
tant  rotational  flow  effects. 

Stability  predictions  for  these  flow  systems  are  usually  obtained  from  mathemat¬ 
ical  models  that  reflect  linear  and  weakly  nonlinear,  inviscid,  irrotational  acoustic 
concepts  (see  Culick,  1990,  Williams, 1985,  Price  and  Flandro,  1993,  Jahnke  and 
Culick,  1993).  The  basic  acoustic  waves  propagate  through  a  fluid  at  rest,  in  the 
first  approximation,  and  are  described  in  terms  of  Fourier  series  of  eigenfunctions 
that  satisfy  '  tally  impermeable  boundary  conditions.  Strictly  speaking  ,  these  so¬ 
lutions  do  not  accomodate  sidewall  injection  nor  exit  plane  flow  at  the  downstream 
end  of  the  cylinder.  As  a  result  the  model  cannot  account  for  an  interaction  between 
the  acoustic  signals  and  a  sidewall  injected  cylinder  flow,  now  known  to  be  a  source 
of  vorticity  generation  (  Flandro  and  Roach, 1992). 

Flandro  (1974)  recognized  that  rotational  flow  effects  play  a  role  in  relatively  thin 
acoustic  boundary  layers  where  viscosity  is  of  significance.  Related  studies  for  inert 
flows  nave  been  carried  out  by  Tien  (1972),  Flandro  (1986),  Hegde  et  al.(1986)  and 
Price  and  Flandro  (1993).  Chemically  and  thermally  active  acoustic  boundary  layer 
flows  are  described  by  Hegde  and  Zinn  (1986),  Sankar  et  al.(1988a,b),  Chen  et  al. 
(1990)  and  Matta  and  Zinn  (1993).  In  these  studies  the  boundary  layer  responds 
passively  to  externally  imposed  disturbances.  The  investigations  are  motivated  by 
a  need  to  understand  how  energy  is  exchanged  between  the  acoustic  disturbances 
and  mean  flow  as  fluid  injected  normally  from  the  wall  is  turned  towards  the  axial 
direction.  Until  recently,  conceptual  understanding  of  this  flow  turning  process  has 
been  based  largely  on  the  viscous  properties  of  the  thin  acoustic  boundary  layer. 

Significant  efforts  have  been  made  to  develop  computational  models  for  acoustic 
boundary  layer  processes.  Baum  and  Levine(1987),  Baum  (1990),  Vuillot  and  Avalon 
(1991)  and  Vuillot  (1991)  have  employed  complete  Navier-Stokes  solvers  to  evaluate 
the  general  cylinder  flow  response  to  imposed  disturbances.  The  last  two  cited  works 
in  particular  demonstrate  that  rotational  disturbances  can  exist  in  a  substantial  por¬ 
tion  of  the  cylinder  for  appropriate  values  of  Re  and  M .  In  other  words,  the  vorticity 
distribution  is  not  always  confined  to  the  traditional,  viscous  acoustic  boundary  layer 
adjacent  to  the  cylinder  wall.  Rather,  for  appropriate  parameter  ranges  it  appears  to 
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be  converted  from  the  wall,  out  into  the  main  portion  of  the  basically  inviscid  cylin¬ 
der  flow  by  the  mean  injected  flow  field.  Clearly,  the  pervasive  presence  of  vortical 
structures  in  the  internal  flow  field  has  significant  consequences  for  the  conceptual 
validity  of  traditional  irrotational  acoustic  stability  models. 

The  presence  of  rotational  inviscid  “acoustic”  disturbances  in  a  laboratory  flow 
field  was  first  implied  by  the  experiments  of  Brown  et  al.  (1986a, b)  and  Brown  and 
Shtieffer  (1992).  The  book  manuscript  of  Price  and  Flandro  (1993)  contains  an  im¬ 
portant  discussion  of  these  issues,  and  several  idealized  mathematical  models  are  used 
to  predict  the  structure  of  basically  inviscid  ,  but  weakly  viscous,  vortical  flow  struc¬ 
tures.  More  recently,  Flandro  and  Roach  (1992)  have  extended  the  theoretical  effort 
to  generate  a  purely  inviscid  model  for  boundary  generated  rotational  disturbances 
that  arise  from  an  inviscid  interaction  between  an  injected  flow  field  and  axial,  planar 
pressure  waves  propagating  in  the  cylinder.  Here  again  the  vortical  structures  are 
convected  out  into  the  cylinder  by  the  radial  component  of  the  injected  velocity  field. 
The  length  scale  of  the  local  vorticity  variation  is  0(M)  smaller  than  the  cylinder 
radius,  implying  that  viscous  shear  stress  may  be  important  on  the  local  level.  This 
work  suggests  that  the  interaction  of  “acoustic  disturbances”  with  a  strongly  injected 
flow  field  leads  to  the  presence  of  important  vortical  structure  in  the  entire  flow  field. 

Zhao  and  Kassoy  (1994)  and  Zhao  et  al.  (1994)  (see  Appendices  (a)  and  (b)) 
provide  an  initial  step  in  formulating  a  rational  mathematical  model  for  internal  flow 
dynamics  which  incorporates  both  acoustic  phenomena  and  vorticity  distributions. 
Perturbation  methods  are  used  to  derive  systematic  approximations  to  the  complete 
compressible  Navier-Stokes  equations.  An  initial-boundary  value  approach  is  used 
to  formulate  a  generalized  unsteady  mathematical  model  capable  of  describing  both 
non-resonant  and  resonant  time  history  of  solutions.  The  boundary  disturbance  is  an 
0{M)  axial,  harmonic  velocity  variation  on  the  closed  end  wall.  The  complete  axial 
velocity  is  found  from  a  superposition  of  three  components  of  equal  magnituide.  First, 
the  steady  component  arises  from  a  solution  to  inviscid,  rotational  Euler  equations 
known  by  Culick  (1966).  Secondly,  there  is  a  planar  irrotational  acoustic  field,  derived 
from  a  traditional  linear  wave  equation  which  satisfies  boundary  conditions  at  the 
closed  and  open  ends  of  the  cylinder.  Finally,  the  rotational,  weakly  nonlinear  viscous 
component  varies  on  two  disparate  length  scales.  The  global  spatial  distribution 
occurs  on  the  radial  length  scale,  while  locally  there  are  velocity  gradients  on  a  scale 
0(M)  smaller.  As  a  result  the  vorticity  of  the  transient  field  is  0{M~^)  larger  than 
that  in  the  steady  field.  The  convection  of  the  vorticity  is  found  to  be  described  by 
a  linear  equation.  Most  significantly,  the  diffusion  of  vorticity  on  the  shorter  length 
scale  is  described  by  a  nonlinear  diffusion  equation.  Solutions  are  found  numerically 
by  using  a  finite  difference  method,  as  well  as  semi-analytically  with  a  Galerkin-like 
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method  developed  by  Wang  and  K«issoy  (1990a,b.c,  1992a, b,  1993). 

Evaluation  and  interpretation  of  the  full  results  show  that  rather  complex  vorticity 
distributions  are  present  throughout  the  cylinder.  Initially,  the  vorticity  is  confined 
to  a  growing  region  between  the  sidewall  and  a  well  defined  front  propagating  toward 
the  cylinder  axis  with  the  local  radial  steady  speed.  The  region  between  the  front  and 
the  axis  contains  the  weakly  vortical  steady  flow  along  with  the  irrotational  acoustic 
solution.  Eventually,  the  entire  chamber  flow  is  filled  with  vorticity  which  vanishes 
only  on  the  axis. 

Fully  computational  methods  are  used  by  Kirkkopru  et  al.(1994)  (see  Appendix  c) 
to  provide  qualitative  supporting  evidence  for  the  solutions  described  in  .Appendices 
a  and  b.  In  this  case  the  driving  disturbance  is  a  harmonic  pressure  transient  applied 
on  the  downstream  exit  plane  of  the  cylinder.  Grid  size  and  spatial  distribution  are 
chosen  to  accomodate  the  multiple  lengthscale  structure  known  from  the  formulation 
exercise.  The  unsteady  rotational  component  of  the  axial  velocity  is  extracted  from 
the  total  value  found  from  a  MacCormack  scheme.  Solution  properties  and  charac¬ 
teristics  are  identical  to  those  found  previously  and  support  the  basic  concepts  of 
vorticity  generation  and  transport. 

Results  obtained  to  date  are  associated  with  single  frequency  harmonic  boundary 
forcing.  Even  then  the  time  response  curves  for  the  total  axial  velocity  are  surprisingly 
complex,  resulting  from  a  combination  of  the  expected  acoustic  response  and  the 
newly  described  rotational  component.  One  may  speculate  that  multiple  frequency 
drivers,  and/or  those  with  greater  complexity  of  spatial  distribution  will  generate 
pointwise  time  response  curves  with  considerable  irregularity,  perhaps  reminiscent  of 
a  turbulent  flow.  In  fact,  it  should  be  recognized  that  the  calculated  or  computed 
rotational  velocity  field  has  certain  characteristics  and  properties  that  are  similar  to 
those  of  a  turbulent  flow. 

One  is  observing  the  evolution  of  a  primarily  inviscid,  rotational,  unsteady  flow 
associated  with  a  large  Reynolds  number.  The  weakly  viscous  effect  has  an  impact 
on  the  shorter  length  scale,  and  the  rotational  flow  distribution  is  weakly  nonlinear. 
Given  that  the  results  are  obtained  from  a  qujisi-analytical  initial  value  analysis, 
the  evolving  flow  field  is  described  much  like  that  in  a  direct  numerical  simulation. 
The  front  separates  a  strongly  rotational  flow  from  one  that  is  characterized  by  an 
irrotational  acoustic  field.  Although  not  shown  in  the  present  work,  we  now  know 
that  if  disturbances  are  generated  directly  by  sidewall  injection  transients  of  sufficient 
spatial  complexity,  the  front  morphology  can  be  quite  corrugated.  As  a  result,  one 
can  find  intermittent  regions  of  strongly  rotational  unsteady  flow  on  an  axial  traverse 
at  a  fixed  radial  distance  from  the  wall.  In  the  present  theoretical  framework  the 
front  configuration  can  actually  be  described  analytically  ! 
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Abstract 

A  mathematicai  model  is  formulated  to  describe  the 
initiation  and  evolution  of  vorticity  in  a  low  Mach 
number  (M),  weakly  viscous  internal  flow  sustained 
by  mass  addition  through  the  side  wall  of  a  long,  nar¬ 
row  cylinder.  An  0(M)  acoustical  disturbance,  gen¬ 
erated  by  a  prescribed  harmonic  transient  endwall  ve¬ 
locity,  interacts  with  the  basically  inviscid  rotational 
steady  injected  flow  to  generate  time  dependent  vor¬ 
ticity  at  the  wall.  The  steady  radial  velocity  com¬ 
ponent  convects  the  vorticity  into  the  chamber.  The 
axial  velocity  associated  with  the  vorticity  fleld  varies 
across  the  cylinder  radius  and  in  particular  has  an 
instantaneous  oscillatory  spatial  distribution  with  a 
characteristic  wave  length  0(M)  smaller  than  the  ra¬ 
dius.  As  the  spatial  structures  are  convected  from  the 
wall  toward  the  axis,  weak  viscous  effects  cause  the 
vorticity  to  diffuse.  The  magnitude  of  the  vorticity 
distribution  in  the  transient  field  is  much  larger  than 
that  in  the  steady  flow  because  the  velocity  gradi¬ 
ents  associated  with  the  former  occur  on  the  smaller 
length  scale. 

An  initial-boundary  value  formulation  is  employed 
and  a  finite  difference  scheme  is  used  to  find  the  non¬ 
linear  unsteady  solutions  when  a  pressure  node  exists 
at  the  downstream  exit  of  the  cylinder.  The  complete 
velocity  consists  of  a  superposition  of  the  steady  flow, 
an  acons^c  (irrotational)  field  and  the  vortical  com¬ 
ponent,  aU  of  the  same  magnitude.  Results  for  the 
velocity  solutions  are  given  when  the  endwall  driving 
frequency  is  nonresonant.  Time  response  output  at 
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a  given  location  shows  that  the  complete  axial  veloc¬ 
ity  field  can  be  quite  irregular,  particularly  if  the  end 
wail  driving  frequency  is  high.  The  formulation  and 
results  provide  a  conceptual  framwork  for  the  study 
of  solid  rocket  engine  chamber  flow  dynamics  which 
supplements  traditional  acoustic  stability  theory  by 
providing  information  about  the  generation  and  evo¬ 
lution  of  vortical  structures. 

1  Introduction 

Predictions  of  flow  d]mamics  in  solid  rocket  engine 
chambers  have  been  obtained  in  the  past  from  acous¬ 
tic  stability  theory,  first  formulated  and  described  by 
Grad^.  Culick  and  co-workers,  (See  Culick^  for  a 
comprehensive  review)  in  particular,  have  developed 
since  mid-1970’s  an  increasingly  sophisticated  linear 
and  weakly  nonlinear  stability  theory  to  describe  dis¬ 
turbance  behavior  observed  in  laboratory  amd  full- 
scale  engines.  Comparisons  of  predicted  pressure 
transients  with  those  in  experiments  appear  to  be 
reasonably  good. 

Recent  experimental  measurements  (Brown  et 
al’~*..  Brown  and  Shaeffer’,)  and  computational  pre¬ 
dictions  (Vuillot  and  Avalon^;  Vuillot^)  of  velocity 
profiles  show  definitive  evidence  of  significant  rota¬ 
tional  flow  structures  distributed  throughout  the  en¬ 
gine  chamber.  In  addition,  Price  and  Flandro*  and 
Flandro  and  Roach*  have  formulated  approximate 
mathematical  models  for  describing  the  generation 
and  evolution  of  vorticity  in  analogies  of  solid  rocket 
engine  chamber.  These  contributions  provide  the  mo¬ 
tivation  for  developing  a  comprehensive  chamber  flow 
dynamics  model  which  includes  both  the  irrotational 
acoustical  phenomena  of  traditional  stability  theory 
and  important  rotational  flow  structure. 

Solid  rocket  engine  stability  predictions  are  usually 
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obtained  from  mathematical  models  that  reflect  lin¬ 
ear  and  'Reakly  nonlinear,  inviscid,  iirotational  acous¬ 
tic  concepts  (Culick*®"^^;  William*®;  Price  and  Flan- 
dro  *).  The  basic  acoustic  waves  propagate  through 
a  fluid  at  rest,  in  the  first  approximation,  and  are 
described  in  terms  of  Fourier  series  of  eigenfunc¬ 
tions  that  satisfy  totally  impermeable  boundary  con¬ 
ditions.  Strictly  speaking,  these  eigenfunctions  do 
not  accomodate  side  wall  injection  (propellant  gasifi¬ 
cation)  nor  chamber  exit  plane  flow  (upstream  of  the 
nossle).  As  a  result  the  model  cannot  account  for  an 
interaction  between  the  acoustic  signals  and  a  side- 
wall  injected  chamber  flow,  now  known  to  be  a  source 
of  vortidty  generation  (Fliuidro  and  Roach*). 

Flandro**  recognised  that  rotational  flow  effects 
play  a  role  in  acoustic  boundary  layers  where  viscosity 
is  of  some  significance.  Related  studies  for  inert  flows 
have  been  carried  out  by  Tien**,  Flandro**,  Hegde 
et  al*^.  and  Price  and  Flandro*.  Chemically  and 
thermally  active  acoustic  boundary  layer  flows  are 
described  by  Hegde  and  Zinn*^,  Sankar  et  al**~**., 
Chen  et  al.**  and  Matta  and  Zinn**.  In  these  studies 
the  boundary  layer  responds  passively  to  externally 
imposed  disturbances.  The  investigation  are  moti¬ 
vated  by  a  need  to  understand  the  energy  exchange 
between  acoustic  disturbances  and  the  mean  flow  as 
the  normally  injected  fluid  is  turned  towards  the  axial 
direction.  Until  recently,  conceptual  understanding 
of  flow  turning  has  been  based  largely  on  the  viscous 
properties  of  the  thin  acoustic  boundary  layer. 

Significant  efforts  have  been  made  to  develop  com¬ 
putational  models  for  acoustic  boundary  processes. 
Baum  and  Levine*®,  Baum*^,  Vuillot  and  Avalon^ 
and  Vuillot*  have  employed  complete  Navier-Stokes 
solvers  to  evaluate  the  general  chamber  flow  response 
to  imposed  disturbances.  The  last  two  cited  works  in 
particular  demonstrate  that  rotational  disturbances 
exist  in  a  substantial  portion  of  the  chamber  flow.  In 
other  words,  vortidty  distributions  are  not  conflned 
to  traditional,  viscous  acoustic  boundary  layers  ad¬ 
jacent  to  the  chamber  wall.  Rather,  they  appear  to 
be  convected  by  the  mean  injected  flow  field  &om  the 
wall  out  into  the  main  portion  of  the  basically  invis¬ 
cid  chamber  flow.  Clearly,  the  pervasive  presence  of 
vortical  structures  in  the  chamber  flow  field  has  sig¬ 
nificant  consequences  for  the  conceptual  validity  of 
traditional  irrotational  uoustic  stability  models. 

The  presence  of  rotational  invisdd  “acoustic”  dis¬ 
turbances  in  a  model  (cold  flow)  chamber  flow  was 
first  inferred  from  the  experiments  of  Brown  et  al.*~* 
and  Sh^mffer  and  Brown*.  The  book  manuscript  of 


Price  and  Fiandro*  contains  am  important  discussion 
of  these  issues,  and  severaJ  mathematical  models  are 
used  to  predict  the  structure  of  basically  inviscid  , 
but  weaddy  viscous,  vortical  flow  structures.  More 
recently,  Flandro  and  Roaich*  have  extended  the  the¬ 
oretical  effort  to  generate  a  purely  inviscid  model 
for  boundary  generated  rotational  disturbances  that 
arise  bom  an  inviscid  interaction  between  an  injected 
flow  field  atnd  axial,  planar  pressure  waves  propagat¬ 
ing  in  the  chamber.  Here  again  the  vortical  structures 
are  convected  out  into  the  chamber  by  the  radial  com¬ 
ponent  of  the  injected  velocity  field.  The  length  scale 
of  the  local  vortidty  variation  is  0(M)  smaller  than 
the  chamber  radius.  This  work  suggests  that  the  in¬ 
teraction  of  “acoustic  disturbances”  with  a  strongly 
injected  flow  field  leads  to  the  presence  of  important 
vortical  structure  in  the  entire  flow  field  of  a  model 
solid  rocket  engine. 

These  rotational  flow  models  also  demonstrate  un¬ 
equivocally  that  the  velodty  field  near  the  wall  is 
a  composite  of  an  axial  planar,  irrotational  acous¬ 
tic  flow  and  a  rotational  flow.  The  short  wave  length 
transverse  spatial  oscillations  of  the  axial  disturbance 
velodty,  seen  in  the  numerical  results  of  Vuillot  and 
Avalon*  as  well  as  more  dearly  in  the  theory  of  Flan- 
dro  and  Roach  *  provide  a  transition  between  the  wall 
and  a  traditional  planar,  axial  acoustic  driving  distur¬ 
bance  present  near  the  axis.  This  suggests  that  the 
concept  of  velodty  coupling  should  be  reconsidered, 
or  perhaps  replaced  by  u  improved  understanding 
of  rotational,  basically  invisdd  disturbance  processes 
adjacent  to  a  burning  propellant  surface. 

The  present  study  is  an  initial  step  in  formulating 
a  mathematical  model  for  solid  rocket  engine  cham¬ 
ber  flow  dynanoics  which  incorporates  both  acous¬ 
tic  phenomena  and  vortical  flow  structures.  A  low 
Mach  number  (M),  large  Reynolds  number  (Re)  in¬ 
ternal  flow,  sustained  by  constant  steady  mass  addi¬ 
tion  through  the  wall  of  a  long  narrow  cylinder,  is 
disturbed  by  an  0(M)  harmonic,  transient  endwall 
velodty.  Acoustic  waves  created  by  the  boundary 
disturbance  interact  with  the  basically  invisdd  rota¬ 
tional  injected  flow  field  to  create  vortidty  on  the 
wall.  The  vortidty  is  convected  into  the  cylinder  by 
the  radial  component  of  the  steady  injected  flow  field 
although  the  flow  structure  is  altered  locally  by  non¬ 
linear  and  viscous  effects. 

Perturbation  methods  are  used  to  derive  system¬ 
atic  approximations  to  the  complete  compressible 
Navier-Stokes  equations.  An  initial-boundary  value 
approach  is  used  to  formulate  a  generalised  unsteady 
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mathematical  model  capable  of  describing  the  com¬ 
plete  time  history  of  solutions.  Fourier  series  repre¬ 
sentations  of  the  velocity  and  pressure  are  obtained  in 
terms  of  eigenfunctions  that  satisfy  all  the  prescribed 
boundary  conditions.  The  complete  axial  velocity  is 
found  from  a  superposition  of  three  components  of 
comparable  magnituide.  First,  the  steady  component 
arises  from  a  solution  to  invisdd,  rotational  Euler 
equations  known  by  Culick^'^.  Secondly,  there  is  a 
planar  irrotational  acoustic  field,  derived  from  a  tra¬ 
ditional  linear  wave  equation.  Finally,  there  is  a  rota¬ 
tional,  weakly  viscous  component  which  varies  on  two 
distinctly  different  length  scales.  The  global  spatial 
distribution  occurs  on  the  radial  length  scale,  while 
locally  there  are  velocity  gradients  on  a  scale  0(M) 
smaller.  As  a  result  the  vortidty  of  the  transient  field 
is  larger  than  that  in  the  steady  field.  The 

convection  of  the  vortidty  is  found  to  be  described 
by  a  linear  equation.  Most  significantly,  the  diffusion 
of  vortidty  on  the  shorter  length  scale  is  described 
by  a  nonlinear  equation.  It  is  demonstrated  that  the 
global  spatial  distribution  on  the  radial  length  scale 
arises  from  the  interaction  between  the  “acoustic  dis¬ 
turbance”  and  the  stron^y  injected  flow  field  present 
in  the  solid  rocket  chamber.  The  nonlinear  and  vis¬ 
cous  effects  are  of  quantitative  importance,but  do  not 
have  much  impact  on  the  overall  flow  structure. 

Solutions  are  found  by  invoking  a  direct  finite  dif> 
ference  scheme  and  a  Galetkin-like  method  used  by 
Wang  and  Kassoy  In  the  latter  cue,  a  cou¬ 

pled,  infinite  system  of  ordinary  differential  equations 
for  Fourier  coefficients  is  truncated  systematically  to 
find  solutions  in  terms  of  a  finite  mode  representa¬ 
tion.  The  two  different  approaches  enable  ns  to  com¬ 
pare  results  for  the  same  cases,  providing  a  means  by 
which  one  can  verify  the  results. 

Evaluation  and  interpretation  of  the  full  results 
show  that  rather  complex  vortical  structures  are 
present  throughtout  the  cylinder  sufficiently  long  af¬ 
ter  the  disturbance  is  initiated  at  the  end  wall.  When 
the  input  from  several  different  Fourier  modes  is  im¬ 
portant,  the  time  response  curves  are  quite  irregular 
and  could  be  mistaken  for  "turbulent”  response. 

2  Mathematical  Formulation 

The  geometrical  configuration  for  a  cylindrical  cham¬ 
ber  with  side  wall  iiyection  is  shown  schematically 
in  Figure  1.  in  which  L'  is  the  axial  length  and  D 
is  the  diameter.  The  radius  R'  is  related  to  D  by 


D  =  2R  .  The  complete  non-dimensional  equations 
describing  the  fluid  dynamics  and  acoustics  for  an 
axis}rmmetric  system  can  be  written  in  the  form 


dt  [r  dr 


)  ^  5(PV.)]  _ 


yM  dr 


ReXdz'^iSi  dz  dr  \ 


1  dp 
yM  dz 


Re  \7  dz"^  dr  \ 

dr^[Bi  dz  dr  \l  Redz'^ 

P  =  pT  (5) 

where 

and  #  is  the  viscous  dissipation  function.  The  non- 
dimensional  variables  are  defined  in' terms  of  dimen¬ 
sional  quantities  (with  a  prime)  by 

p  =  i,P  =  ^,T=^,Vr  =  ^, 


=  jr,r  - 
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k  =  =  (6) 

The  reference  value  pQ  measures  the  basic  static  pres¬ 
sure  in  the  chamber,  while  the  analogous  density 
and  temperature  values  Po.  ^0.  respectively,  repre¬ 
sent  properties  of  the  fluid  injected  from  the  walls. 
The  characteristic  injection  speed  V^q  is  related  to 
the  characteristic  axial  speed  V^q  by  the  approximate 
mass  conservation  relationship  =  6V^q  where  the 

aspect  ratio  6  =  Characteristic  length  scales  for 
the  axied  and  radial  dimensions  are  defined  by  L'  and 
R  respectively.  Time  is  nondimensionalized  with 
respect  to  the  sixial  acoustic  time  t^  =  ^,  where 

C'q  =  (7^  Tg)^  is  the  characteristic  sound  speed. 
The  reference  material  properties  ibg,  /ig  and  ^*0  are 
defined  at  temperature  Tg.  The  parameter  7  is  the 
ratio  of  specific  heats  and 


Mo  *0  ^0 


(7) 


where  typically  the  Prandtl  number  Pr=  0(1),  the 
axial  Mach  number  M  C  1  and  axial  Reynolds  num¬ 
ber  Re  >  1  in  the  chamber. 

Initially,  the  steady  flow  in  the  system  is  driven  by 
a  spatially  distributed  normal  injection  from  the  wall 
where  the  no-slip  condition  is  satisfied.  Symmetry 
prevails  along  the  axis  at  the  cylinder.  The  closed 
end  of  the  cylinder  is  impermeable,  and  the  exit  is 
assumed  to  be  a  pressure  node.  The  latter  assump¬ 
tion  is  a  mathematical  modelling  convenience.  The 
boundary  conditions  may  be  written  as 

r  =  0;  V;  =  ^  =  0, 
r  =  1;  K  =  -Vr^{z),  V,  =  0,  r  =  1  (8) 

X  =  0;  V.  =  0,  2  =  1;  P  =  1.  (9) 


The  injection  rate  is  characteristic  of  the  gasifi¬ 
cation  velocity  created  by  solid  propellant  burning. 

Large  injection  prevails,  in  the  sense  that  >  ii,. 

-  Ret 

(Cole  and  Aroesty^®),  which  implies  that  Re  »  6^, 
where  5  1  is  used  to  describe  a  long,  narrow  cham¬ 

ber. 


3  Steady  State  Flow 

The  steady  state  flow  generated  by  time  independent 
mass  addition  on  the  side  wall  can  be  described  in 
terms  of  the  variables: 


(P)  P.  T)  ~  1  Af*(Po»,  Po».  Tg,)  , 

(V..K)~(V'.0.,V;g.)+..-  (11) 


Equation  (11)  can  be  used  in  (l) — (6)  to  find  the 
leading  order  equations  valid  for  the  limit  M— >0,  de¬ 
scribing  an  inviscid  rotational  flow  that  satisfies 


VrOt  +  v,o,  - 


1  8P0, 
7  dz  ' 


(12) 

(13) 

(14) 


and  boundary  conditions  in  (8)-(9).  The  transport 
terms  ue  excluded  from  the  leading  order  equations 
because  ^  >  1.  Equation  (13)  arises  because  the 
aspect  ratio  ^  >  1.  The  solutions  for  the  radial  and 
axial  velocity,  as  well  as  the  pressure  distribution  can 
be  written  in  the  form 


Vro.  =  -^un{Zr^), 

Vmo,  =  ("-/o  V'rw(r)dr)  cos(f  r*),  (15) 


Plow  disturbances  are  created  at  z  =  0  by  imposing 
a  harmonic  end  wall  disturbance  in  the  axial  velocity 
which  is  independent  of  the  radial  coordinate; 

z  =  0;  V,  =  Asinu/t,  t  >  0;  0  <  r  <  1,  (10) 

where  the  amplitude  A=0(1). 

It  should  be  noted  that  the  imposed  disturbance  on 
the  left  end  wall  is  comparable  in  magnitude  to  that  of 
the  steady  axial  speed.  As  a  result  this  fluid  system 
will  be  weakly  nonlinear,  and  the  theory  developed 
here  differs  from  traditional  linear  small  disturbance 
theory. 


Po.  =  Si  [v;.(i)  So  VrAr)dr\  dz,  (16) 

where  — 14«(z)  is  an  arbitary  time-independent  side 
wall  injection  distribution.  Related  solutions  can  be 
found  in  Culick^°,  Price  and  Flandro*  and  Flrmdro 
and  Roach^. 
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4  Rotational  Flow  Across  the 
Cylinder 

4.1  Problem  statement 

The  objective  now  is  to  develop  &  mathematical 
model  fot  the  disturbance  evolution  when  weak  rota* 
tional  effects  are  significant  across  the  entire  cylinder. 
The  asymptotic  expansions  for  the  velocity  compo¬ 
nents  and  thermodynamic  variables  are 

V.  ~  Ko.(*,r)  +  5;  (17) 

n=0 

K~Ko.(s,r)-^53if»t;„(z,r.t)  (18) 

nsO 


Af  — •  0.  First,  the  spatially  homogeneous  boundary 
forcing  in  (10)  and  the  condition  ^  >  1  imply  that 
VfQ  =  0.  Then, 


dRo  dRo 
at  5ra 

dV,o  _ 

(23) 

dv,o  _  laPo 

at  dr-i  y  dz 

(24) 

8P0  _  aPo  _  Q 

dvi  dvi 

(25) 

aso  a^Q  _  (7  ~  1)  dPo 

dt  dr2  7  9t 

(26) 

Po  =  Rq  +  So- 

(27) 

(p,p.r)~i  +  M5;j;f*‘(p„,p,.«,)  (19) 

n=0 

It  is  recognised  that  two  disparate  length  scales  are 
important;  the  tube  radius  and  a  much  shorter  length 
associated  with  the  radial  distance  traveled  by  a  fluid 
particle  on  the  acoustic  timescale.  A  multiple  scale 
analysis  will  be  carried  out  in  terms  of  the  variables 
rj  and  rj  defined  by 

..  =  1-,,  (20) 

The  second  transformation  includes  an  integral  of  the 
steady  radial  velocity  field  for  the  case  of  constant 
steady  wall  iiyeetion  s  1.  The  integral  trans¬ 
formation  simplifies  the  describing  equations  consid¬ 
erably.  It  is  noted  that  when  the  center  line  is  ap¬ 
proached,  1,  the  integral  diverges  and  rj  — »  00. 

The  partial  derivatives  with  respect  to  r  in  equa¬ 
tions  (l)-(5)  most  be  replaced  by 


(^)  -  {hi)  -  {mL.  ar"er,} 


4.2  Lowest  order  mathematical  model 

The  relations  (7) — (19)  can  be  substituted  into  (1) — 
(5)  to  find  the  leading  order  equations  in  the  limit 


FoUowing  a  procedure  related  to  that  described  by 
Lagerstrom^,  and  similar  to  that  employed  by  Price 
and  Flandro*  and  Flandro  and  Roach*,  the  variables, 
except  for  Pqi  ue  divided  into  irrotational  planar  and 
rotational  nonplanar  parts, 

VMO^Vl'of(z,t)-hWo(z,t,ri,r2),  (28) 

Jlo  =  Rof(z,t) -(- Ro(^.t,ri,ra),  (29) 

=  ^ov(*it) +  ra).  (30) 

When  (28)— {30)  are  used  in  (23) — (27),  the  planar 
equations  are  found  to  be  in  the  form 

at  ~  dz  ' 

(32) 

at  ~  y  dz'  ^  ’ 

96op  _  7  -  1  9Po 

at  y  at  '  ^  ^ 

Po  =  Ror  +  ^Op  (34) 

The  initial/bonndary  conditions  are  given  by: 

t  =  0,  Wo,  =  0,  ^=0.  (35) 

i  =  0,  Wqp  =  sin  (i/t ;  (36) 


*  =  1,  ^  =  0-  (37) 

The  nonresonant  acoustic  solution  for  the  planar 
contribution  is 

Wop{t, z)  =  ^(u>t) 
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-  siii(A„f))sin(A„2),  A„  w,  (38) 


Po{t,i)  =  ^a„(-cos(u;t) 

n  =  l 

+  cos(A„t))cos(A*2),  A„  ^  w,(39) 

^0  =  yRop  (40) 

where  o*  =  -  and  =  (n  -  A)t.  The  first 

terms  in  the  sums  of  (38)  nnd  (39)  arise  from  the  forc¬ 
ing  at  frequency  w,  and  the  second  terms  describes 
the  eigenfunction  response.  Only  the  nonresonant 
case  will  be  considered  in  the  present  work. 

The  equations  for  the  rotational  components  are 

at  art  dz  ~  \VroJ  drj  ’ 

(«> 

^  +  Sq  s  Q.  (44) 

Equations  (43)  and  (44)  can  be  combined  to  show 
that  the  leading  order  vortical  flow  is  incompressible: 

^  +  («) 

Therefore,  equation  (41) — (44)  can  be  rewritten  in 
simpler  forms 


*  VKoJ  an 


awo  awo 


Ro  +  So  —  (49) 

and  must  be  solved  with  respect  to  the  initial  and 
boundary  conditions: 


t  =  0;  Wo  =  Q, 


2  =  0;  Wa  -•  0, 


^=0, 

oz 


,  dWo  ^  dOo 

r,  =  ,.  ^0.  _=0,  (53) 


ri  =  rj  =  0;  Wq  =  -f'"  p(t,  2),«o  =  -6op(t,  z)  (54) 

Equation  (54)  corresponds  to  the  no-sUp  condition 
and  isothermal  flow  injection.  Equations  (42),  (43) 
and  (45)  show  that  Wq,  Sq  and  ^  are  invariaint  cn  a 
characteristic  line  defined  by 

n^t-n,  (55) 

but  vary  across  the  q  lines.  The  r,  =  0  line  enters 
the  chamber  at  t  =  0  through  side  wall  (rj  =  0)  and 
subsequently,  at  t  =  c,  q  =  c  appears  at  n  =  0.  At  a 
particular  time  T,  constant  q  lines,  which  range  from 
0  to  T,  are  transported  toward  the  *Ti«  by  convection 
at  the  local  radial  steady  velocity. 

The  inviscid  equation  in  (42)  can  be  combined  with 
the  first  of  (54)  to  show  that  the  wall  is  a  source  of 
vorticity  because 

9Wo(t,2,0,0)  _  awo  _  awo,(t,z) 

arj  ~~  at 

where  Wo,  is  known  &om  (38).  It  is  noted  that 
the  largest  unsteady  nondimensional  vorticity  term 
is  given  by  O*  =  (jp^)^.  The  parameter, 
arises  &om  large  gr^ents  occuring  in  the  spatially 
oscillatory  velocity  profile  on  the  short  length  scale 
n-  Equation  (47)  also  shows  that  the  vorticity  gen¬ 
erated  at  the  wall  is  convected  into  the  cylinder  by 
the  steady  radial  velocity  field  V'rOt(''). 

In  general  solutions  to  the  inviscid  first  order  hy¬ 
perbolic  equations  in  (47)-(49)  can  be  written  sym¬ 
bolically  u. 


Wo  =  Wo[r),n,z), 


7)  =  t-n 


ao  =  So(v,  >’1,  z)  =  -Ro,  (58) 

In  order  to  derive  the  explicit  functional  dependence 
of  the  variables,  higher  order  equations  need  to  be  in¬ 
vestigated.  Once  I^o  is  found,  then  the  mass  conser 
vation  equation  (46)  can  be  integrated  with  respect 
to  r2  to  find  the  radial  velocity  V^i.  The  vortical  tem¬ 
perature  and  density  fields  can  be  found  using  related 
methods. 
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4.3  Higher  order  consideration 

Equations  (17) — (19)  can  be  combined  with  (1) — (5) 
to  find  the  0(M)  equation  set.  The  same  procedure 
used  to  find  the  leading  order  solution  is  employed 
so  that  the  variables,  except  for  Pi,  are  divided  into 
irrotational  planar  and  rotational  nonplanar  parts, 

V,i  =  Wi,{z,t)  +  Wi(z,t,ruT2),  (59) 

Ri  =  +  Ri{i,i,ri,r2),  (60) 

=tfi,(z,t)  +  fii(s,t,ri,rj).  (61) 

The  planar,  acoustic  equations, 

dRi,  ^  djRopWop)  ^ 


If  the  transformation  of  the  coordinate  system  from 
(t,»'i.'‘2.-t)  to  in,  ri,r2,  z)  is  made,  then  the  derivatives 
with  respect  to  t  and  r2  must  be  replaced  by 


It  follows  that  (66)  can  be  written  as 

dWi  _  1  d^_  djV.o,  +  Wq) 

^’■2  ^ 


1  ^(Kliyo)  _  djVrlW..) 

VrOi  dri  d-rj 


1  3(P,  -  Pa,) 
7  9z 


1  _  dPo  .  .. 


-Poi7  -  1)^.  (64) 

Pi  =  Rip  +  9ip  +  Rop^Op  (65) 

are  not  considered  further  here. 

The  largest  possible  viscous  effect  occurs  when 
6^/Re=M^  and  leads  to  the  higher  order  rotation^ 
equation  for  axial  momentum 


dWi  dWi 
dt  ^  drj 


1  1  f.  dPo 


-(V.O,  +  V,„) 


a(v'.o, +r,„) 


+VrO. 


d{V,o.  +  V,o) 


VVroJ  dr  2 


'Vro.^.  (66) 

ari 

It  is  (66)  that  will  give  us  information  about  the 
behaviour  of  the  leading  order  axial  speed  solution. 


OTi  7  dz 

An  integration  of  (69)  with  respect  to  rj,  holding  r), 
ri  and  z  fixed  will  generate  secular  growth  in  r;  un¬ 
less  certain  terms  are  suppressed.  In  considering  the 
impact  of  each  term,  it  is  important  to  remember  that 
the  harmonic  t>dependence  of  the  planar  acoustic  so¬ 
lutions  in  (38)>(40)  must  be  rewritten  in  terms  of  ij 
and  rj  by  using  the  transformation  r;  =  t  -  rj  in  (55). 
When  written  in  the  coordinate  system  (z,t,  ri.rj) 
the  suppressed  terms  take  the  form; 

I  d^Wo  ,  ,,  dWo  ,„i,  dWo 

-lVo^^-K,o.^=0,  (70) 

which  is  a  nonlinear  diffusion  equation  for  the  rota¬ 
tional  arial  velocity  Wq  with  a  time-like  variable  ri. 
The  solution  for  Wq  must  satisfy  an ‘initial”  condi¬ 
tion; 

W^o(f,  2!.  »'i  =  0.  rz)  =  -V^opin,  *)  »7  >  0 

=  ^  En=l  “»  «“(‘"7?) 

-  sin(A„»7))  sin(A«z),  A„  ^  w, 

Wo(f.2,»-i  =  0,r2)  =  0  f7<0  (71) 

and  a  boundary  condition  at  the  center  line: 

rj  — ►  oo,  =  0.  (72) 
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In  addition  a  condition  must  be  specified  on  > 
0,r2  =  0  which  is  compatable  with  (71)  at  the  point 
ri  =  r2  =  0.  This  is  necessary  because  ri  and  r2 
are  treated  as  independent  viiriables  in  (70).  The 
reasonable  choice  is  given  by 

Woit,z,rur2  =  0)= -Wopit,z).  (73) 


where 


Onn,n,  =  2A.„,  /  sin(A„, z) cos(A„,2) sin(A„2)d2  (79) 
Jo 

^nni  “  Afij  /  2Cos(A„,2)sin(An2)d2  (80) 

Jo 


The  nonlinear  term  in  (70),  WqWo,  is  present  in  our 
problem  because  the  0(M)  boundary  disturbance 
is  Ivger  than  that  used  in  earlier,  basically  linear 
studies  (Flandro  iuid  Roach^).  Its  presence  suggests 
that  wave  steepening,  acoustic  streaming  and  other 
forms  of  instability  may  occur  in  the  evolving  flow 
held.  If  the  imposed  piston  driven  disturbance  is 
smaller,and/or  axial  variations  are  ignored,  then  a 
linear,  viscous  diffusion  equation  is  derived,  which  is 
related  to  that  used  by  Price  and  Flandro*. 

Given  the  forcing  condition  in  (71),  one  approach 
to  finding  the  solution  to  (70)  is  based  on  using  the 
eigenfunction  set  {sin(An2)},n  =  1, 2,  •  •  •; 

00 

Woit,  z,  ri,  rz)  =  ^  Anit,  ri,  r2)  sin(A„2)  (74) 

.  nsl 

Coupled  partial  differential  equations  for  the 
Fourier  coefficients  Ant  ue  found  by  using  (74)  in 
(70)  and  invoking  orthogonality  conditions  for  the 
eigenfunction  set  on  the  interval  [0,1].  The  results 
are; 

1  5* An  ,  aAn  1  dV,o,  ,  2  /  ^ 

drl  ^  dn  Vro.  dz  ^  Vro. 

y  <liinjna.4n,  An,  +  ( — hnnjAn,  j  ~  ® 

and  are  subject  to; 

(a)  Initial  Condition: 

An(t,»-i  =  0,r2)  =  ^On  «nw(^  -  ^2) 

-  sin  A„(t  -  rz))  0  <  rz  <  t 

=  0  r2  >  t  (76) 

(b)  Boundary  Conditions: 

’•2  — 00,  ^=0(77) 
An(t,ri,r2  =  0)  =  ^On  Sin(ci;t)  -  sin(A„t)^  (78) 


Equation  (75)  is  a  coupled  system  of  quasi-linear  dif¬ 
fusion  equations  with  a  time  like  variable  rj ,  and  non¬ 
linear  source  terms.  The  physical  time  t  is  a  parame¬ 
ter  of  the  differential  equations,  appearing  explicitly 
only  in  the  boundary  conditions.  Solutions  to  (75)- 
(80)  are  described  in  Section  6,  based  on  a  system 
truncation  approach  to  find  a  finite  number  of  An’s. 

A  second  solution  approach  is  based  on  a  fi¬ 
nite  difference  approximation  to  (70)-(73).  Solu¬ 
tions  for  Wo(2,  t,  ri,  rz)  are  found  by  employing  a  sec¬ 
ond  order  accurate  Adam-Bashforth/Crank-Nicolson 
scheme.  It  is  relatively  easy  to  implement  and  can  be 
used  to  test  the  accuracy  of  the  truncated  mode  so¬ 
lution.  Most  of  the  results  presented  here  have  been 
found  in  this  way. 

5  Finite  Difference  Solutions 
for  the  Nonlinear  System 

Solutions  to  (70)-(73)  for  Wo{z,t,ri,r2)  have  been 
found  by  using  a  finite  difference  method  based 
on  the  second  order  accurate  0(5rf,6r^)  Adam- 
Bashforth/Crank-Nicolson  scheme.  The  real  solution 
lies  on  the  locus  defined  by  (20)  relating  ri  and  rz. 
In  general  this  curve  does  not  coincide  with  the  com¬ 
putational  grid  points.  Hence  cubic  spline  interpola¬ 
tions  in  rz  are  employed  at  every  chosen  ri  value  to 
obtain  the  desired  real  solution.  The  far  end  bound¬ 
ary  condition  for  rz  — •  00  is  implemented  by  pro¬ 
viding  a  sufficiently  large  number  of  grid  points  be¬ 
tween  the  convected  outer  edge  of  the  rotational  layer 
and  the  finite  location  of  the  computational  bound¬ 
ary  with  respect  to  rz.  A  total  of  1000  grid  points 
in  the  rz  direction  with  £r2=0.1  is  chosen.  The  func¬ 
tion  and  derivatives  must  remain  sero  at  a  significant 
number  of  nodes  in  order  to  ensure -that  conditions 
at  the  computational  boundary  do  not  constrain  the 
solution. 

At  each  value  of  the  “parameter”  t,  the  integration 
is  initiated  with  the  initial  conditions  in  (71),  subject 
to  the  boundary  conditions  in  (72)  and  (73)  .  The 
spatial  distribution  of  the  solution  with  respect  to  rz 
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evolves  as  the  “time'like”  variable  ri  increases.  Inte¬ 
gration  is  carried  out  to  a  sufficiently  large  value  of  ri 
to  ensure  that  adequate  data  fields  Wq^z,  t,  r^,  rj)  are 
available  on  the  locus  curve  relating  rx  and  rj  defined 
by  (20).  Then  the  actual  solution  is  found  from  the 
intersection  of  the  surface  defined  by  IVo(z,t,  ri,r2) 
and  the  vertical  plane  from  the  locus  curve  relating 
ri  and  rj. 

The  first  test  case  studied  is  for  «>  =  1.0,  which  is  a 
nonresonant  &equency  smaller  than  the  first  natural 
frequency  Ax  =  |.  This  particular  case  allows  ns 
to  develop  a  relatively  simple  solution  with  Tniwimal 
computational  time. 

It  is  important  to  resolve  the  solution  in  the  axial 
direction  by  choosing  a  sufficiently  large  number  of 
grid  points  in  the  s-direction,  Solution  com¬ 
parisons  for  Km»x  =  21  and  37  at  t  =  40,  MsO.Ol 

and  Z=0.2S,  0.5  and  0.75  show  considerable  differ¬ 
ence  for  the  smallest  number.  Excellent  comparisons 
for  =  21  and  37  suggest  that  the  former  is  ad¬ 
equate  for  fa;  =  1.0  . 

Figures  2,3  and  4b  give  results  for  the  rotational  ax¬ 
ial  velocity  Wq  as  a  function  of  rx  at  t=0.5  when  t= 
20,30  and  40.  At  each  time,  the  solution  consists  of 
several  spatial  oscillations  and  a  thin  repon  of  expo¬ 
nential  damping  near  a  specific  radial  location  rie(t) 
beyond  which  vortidty  is  exponentially  small.  For  a 
given  value  of  t,  one  first  finds  rje  from  17  =  0  s  t^rze 
and  then  uses  the  inverse  of  (20)  to  determine  rxe.  It 
should  also  be  noted  &om  the  definition  of  17  in  (55) 
and  the  second  of  (20)  that  the  front  moves  at  the 
nondimensional  speed  =  —MV,o,{r)  which 

corresponds  to  the  steady  radial  speed  in  dimensional 
terms.  The  &ont  speed  vanishes  as  the  cylinder  axis 
is  approached.  At  t  =  40  the  rotational  flow  distri¬ 
bution  in  Figure  4b  extends  out  about  44%  of  the 
cylinder  radius. 

The  corresponding  scaled  vortidty  distribution, 
calculated  from  ^  given  in  Figure  5.  The 

dimensionless  vortidty  0#,  defined  below  (56)  is 
0{M~^).  One  should  note  the  significant  magnitude 
of  the  spatial  variations  in  vortidty  given  the  scaling 
factor  ■  The  large  amplitude  is  assodated  with 
the  0(MjTength  scale  of  the  axial  speed  gradient  and 
the  factor  3f~‘. 

The  spatial  distribution  of  the  axial  rotational  ve¬ 
locity  at  each  time  can  be  explained  in  physical  terms 
by  considering  the  interaction  between  the  propar 
gating  planar  acoustic  wave  solution,  described  by 
(38)  and  (39),  and  the  steady  injected  flow  field  in 


(15).  Fluid  partides  injected  normally  from  the  wall 
with  K*  =  1  at  a  si>ecified  z-location  experience  an 
approximately  harmonic  variation  in  the  local  pres¬ 
sure  gradient  given  by  (39).  During  periods  of  neg¬ 
ative  (positive)gradients  the  partides  are  accderated 
downstream  (upstream).  As  a  result,  near  the  wall 
one  will  observe  alternating  periods  of  positive  and 
negative  axial  vdodty.  The  steady  raffial  velodty 
field  carries  these  alternating  regions  of  forward  and 
reverse  flow  away  from  the  wall  toward  the  axis.  Part 
of  the  fluid  partide  response  is  purdy  acoustic.  The 
remainder  is  given  by  Wq.  In  this  sense  the  spatial 
pattern  of  the  rotational  axial  vdodty  at  fixed  t  as 
shown  in  Figures  2,3  and  4b  reflects  the  historical 
behaviour  of  the  local  pressure  gradient  at  the  *ria] 
location.  The  0(M)  length  scale  of  the  transverse 
spatial  oscillatioiu  can  now  be  explained  easily  since 
the  approximatdy  harmonic  pressure  gradient  varia¬ 
tion  occurs  on  the  acoustic  time  scale  (for  fa;=0(l)) 
during  which  only  limited  radial  motion  is  possible. 

Two  lulditional  results  should  be  noted.  The  wave 
length  of  the  spatial  oscillations  decrease  as  rx  in¬ 
creases  toward  the  axis.  This  occurs  becanse  the 
v(^rtidty  front  speed  (the  steady  radial  vdodty)  de¬ 
clines  as  rx  increases.  Finally,  a  comparison  of  Fig¬ 
ures  4a,b,c  demonstrates  that  the  rotaional  axial  ve¬ 
lodty  dependes  on  the  axial  location  as  weU  as  the 
radial  location. 

Figures  6  shows  a  *ffinear"  solution  for  the  rota¬ 
tional  axial  vdodty  *  function  of  rx  at  t=40, 

M=0.01  and  s=0.5  obtained  by  solving  (70)  with 
nonlmear  term  reduced  by  a  factor  of  10~^. 

Jfmss  =  21  is  used  in  the  calculation.  This  calcu¬ 
lation  is  done  to  assess  the  impact  of  the  nonlinear 
convection  term  on  the  spatial  structure  of  the  so¬ 
lution  in  Figure  4b.  A  comparison  shows  that  the 
nonlinear  term  has  a  quantitative  effect  on  the  spa¬ 
tial  distribution  of  Wo,  but  does  not  fundamentally 
control  the  qualitative  spatial  oscillations. 

Figure  7  is  the  counterpart  to  Figure  4b  with  a 
reduced  viscous  effect.  In  this  case  (70)  is  solved  with 
the  vsKons  term,  multiplied  by  a  factor  of 

0.5.  The  basic  patterns  of  the  complete  solution  in 
Figure  4  persist.  There  exist  some  differences  in  the 
amplitudes  of  the  curves  in  these  figures  that  can  be 
attributed  to  the  viscosity  reduction. 

Figure  8  and  Figure  9  are  the  counterparts  of  the 
nonlinear  result  in  Figure  4b  with  M=0.05  and  0.1  re¬ 
spectively,  and  t=40.  These  Mach  numbers  are  more 
reflective  of  the  flow  conditions  in  a  rocket  engine. 
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One  out  see  that  the  larger  Mach  number  (stronger 
blowing)  is  associated  with  a  thicker  unsteady  vorti¬ 
cal  region.  The  vorticity  has  filled  the  entire  chamber 
for  M=0.05  and  for  M=0.1  at  t=40  while  only  44%  of 
the  chamber  is  filled  in  Figure  4  with  M=0.01.  This  is 
expected  because  of  faster  radial  convection  of  vortic¬ 
ity  by  the  steady  radial  velocity.  It  should  be  noted 
that  the  local  velocity  gradients  are  smaller  in  Fig¬ 
ures  8  and  9  so  that  the  magnitude  of  the  unsteady 
vorticity  is  similarly  reduced  for  higher  Mach  number 
systems.  Larger  velocity  gradients  can  be  obtained  at 
a  given  Mach  number  for  larger  forcing  frequency  u. 

Figure  10  provides  results  at  ri  =  0.2  for  the  time- 
variation  of  the  planar  acoustic  axial  speed  Wijp^  the 
rotational  axial  speed  Wq  and  their  sum  defined 
in  (28),  at  s=0.5  with  Kn»x  —  31  and  M=0.01.  It 
can  be  seen  from  Figure  10b  that  the  curve  for  the 
rotational  solution  resembles  the  planu  acoustic  re¬ 
sponse  in  Figure  10a  but  differs  in  amplitude  and 
phase.  At  ri  =  0.2  the  rotational  response  appears 
after  a  delay  of  almost  18  axial  acoustic  time  units, 
the  time  needed  for  the  vorticity  wave  front  initiated 
at  the  wall  to  travel  out  to  the  specified  radial  loca¬ 
tion.  It  is  noted  that  the  amplitude  of  the  rotational 
response  is  a  little  smaller  than  that  of  the  acoustic 
solution.  At  the  location  ri  =  0.2,  phase  differences 
are  relatively  small  and  the  sum  in  Figure  10c  shows 
a  total  response  of  significant  amplitude.  This  am¬ 
plitude  will  actually  increase  as  ri  decreases  until  lo¬ 
cations  very  close  to  the  wall  are  reached,  then  the 
no-slip  conditions  prevails  as  — »  0. 

The  smaller  amplitude  vortical  radial  velocity  can 
be  calculated  by  integrating  equation  (46).  Figure 
11  describes  the  rotational  radial  velocity  Vri  as  a 
function  of  ri  at  t=40,  M=0.01  and  s=  0.5  .  Beyond 
the  vorticity  front  the  transient  radial  velocity  dis¬ 
tribution  is  a  known  constant  multiplying  the  steady 
radial  velocity. 

The  second  case  studied  is  for  a  driving  frequency 
a;  =  2.5  .  The  solution  for  Wq  at  t=40,  M=0.01  and 
z=0.5  is  given  in  Figure  12  for  Km*x  —  21  .  The 
absolute  maximum  amplitude  of  the  w  =  1  solution 
at  s=0.25  is  about  4.  In  comparison,  the  analogous 
queuitity  is  a  quarter  as  large  for  the  case  u  =  2.5. 

The  graphs  show  that  the  spatial  distribution  and 
the  time- response  curves  for  <<;  =  1.0  and  2.5  have 
similar  characteristics.  The  curves  for  w  =  2.5  are 
more  irregular  than  their  u  =  1.0  counterparts  be¬ 
cause  the  driving  frequency  is  higher.  The  relative 
complexity  for  u  =  2.5  suggests  that  if  a  multidimen¬ 
sional  irrotational  acoustic  field  were  present,  involv¬ 


ing  many  Fourier  modes,  the  vortical  time-response 
might  be  quite  irregular.  In  this  sense,  one  could 
ask  whether  ‘Hurbulent”  responses  observed  in  solid 
rocket  chamber  models  are  caused  by  wall  generated 
vorticity  that  is  convected  into  the  chamber  by  the 
injected  field. 

6  Modal  Solutions  to  the  Non¬ 
linear  System 

Solutions  found  from  (74)-(80)  enable  one  to  demon¬ 
strate  how  the  energy  is  distributed  among  different 
Fourier  modes.  The  coefficients  An(ei,r2),  n=l-N 
for  specified  N,  have  been  found  by  using  a  finite  dif¬ 
ference  method  based  on  the  same  second  order  ac¬ 
curate  0(6r\y6T\)  Adam-Bashforth/Crank-Nicolson 
scheme.  The  same  procedure  applied  in  finding  the 
solution  of  Wq  on  the  real  locus  is  applied  to  each  one 
of  the  An(*‘i,  r:)  to  find  the  solution  of  A^s  on  the 
real  locus.  Once  the  A'^s  are  known,  the  solution  for 
Wq  is  found  by  summing  up  N  modal  contributions 
according  to  (74).  Careful  attention  must  be  given  to 
the  value  of  N  in  order  to  assure  that  the  solutions 
are  sufficiently  accurate.  In  particular,  solutions  for 
W'o  based  on  N  and  N  +  L  modes,  L  >  0,  are  com¬ 
pared  until  the  results  exhibit  no  appreciable  change 
for  an  incremental  L  value.  In  Figures  13  and  14,  6 
and  10  mode  solutions  for  >*1  M=0.01,s=0.5 

and  t=40  for  u  =  1.0  have  been  compared  with  each 
other  to  ascertain  that  the  former  are  sufficiently  re¬ 
solved,  implying  that  energy  are  concentrated  in  the 
first  6  modes.  The  comparison  of  Figures  13  and 
14  with  Figure  9b  reveals  that  there  is  little  differ¬ 
ence  between  the  spatial  distribution  curve  from  the 
6  mode  solution  and  that  from  the  direct  finite  differ¬ 
ence  solution  of  nonlinear  equation  with  Kmax  =  21. 
This  is  no  coincidence  because  the  grid  sise  in  the  ax¬ 
ial  direction  for  Ktnax—^^  fr  which  is  just  small 
enough  to  resolve  the  first  six  Fourier  modes  in  ax¬ 
ial  direction  according  to  {Az)maz  =  ^here 

^  =  (»*  -  ?)»• 

The  satisfying  comparison  between  the  finite  differ¬ 
ence  and  mod^  solutions  provide  strong  confidence 
in  the  characteristic  solution  properties.  Solutions  to 
the  former  system  are  less  expensive  in  terms  of  CPU 
time  and  more  complete  if  a  sufficient  number  of  grid 
points  are  provided  in  s  direction  .  Solutions  to  the 
latter  system  are  able  to  demonstrate  the  energy  dis¬ 
tribution  among  different  modes. 
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7  Conclusions 

Systematic  mathematical  methods  have  been  em¬ 
ployed  to  foimnlate  a  mathematical  model  for  the 
creation  and  evolution  of  rotational  flow  in  an  ideal¬ 
ised  solid  rocket  engine  chamber.  Boundary  driven 
axial,  planar  acoustic  waves  interact  with  an  invis- 
cid,  rotational ,  injection  induced  steady  flow  to  pro¬ 
duce  time  dependent  vortidty  at  the  sidewall  of  the 
cylinder.  The  vortidty  is  convected  into  the  entire 
chamber  by  the  steady  radial  velodty  fleld. 

In  contrast  to  traditional  acoustic  stability  theory, 
solutions  are  developed  in  the  context  of  an  initial- 
boundary  value  problem,  in  order  to  study  the  nat¬ 
ural  evolution  of  the  flow  disturbances  arising  &om 
specified  boundary  fordng.  This  study  parallels  typi¬ 
cal  computational  studies  of  similar  problems  (Baum 
and  Levine^*;  Vuillot  and  Avalon^;  Vuillot^).  The 
solutions  found  in  this  work  describe  nonresonant 
boundary  driven  disturbances  on  the  axial  acoustic 
time  scale  of  the  chamber. 

The  analysis  demonstrates  that  the  chamber  flow  is 
basically  invisdd  and  rotational.  The  axial  velodty 
gradients  assodated  with  radial  spatial  oscillations 
occur  on  a  short  length  scale,  0(M)  less  than  the 
radius.  Hence  the  unsteady  vortidty  distribution  in 
the  chamber  is  the  dominant  rotational  flow  effect 
in  comparison  to  the  steady  invisdd  rotational  flow 
caused  by  steady  side  wall  iigection.  The  basic  un¬ 
steady  vortidty  structure  arises  from  an  interaction 
between  the  "acoustic  disturbance”  and  the  strongly 
injected  flow  field.  The  nonlinear  and  viscous  effects 
on  the  0(M)  length  scale,  affect  the  specific  veloc¬ 
ity  distribution,  but  have  little  impact  on  the  overall 
structure. 

Solutions  are  sought  from  the  duect  finite  differ¬ 
ence  calculation  of  equation  (70)  and  of  the  Fourier 
mode  presentation  (75)  compatible  with  a  down¬ 
stream  pressure  node  boundary  condition.  The  for¬ 
mer  calculation  requites  suifident  number  of  grid 
points  in  the  axial  direction  in  order  to  resolve  suf¬ 
ficient  number  of  spatial  mode  in  the  axid  direc¬ 
tion.  The  second  system  is  truncated  to  find  a  lim¬ 
ited  number  of  Fourier  modal  amplitudes.  As  long  as 
the  transient  flow  system  is  nonresonant  (the  fordng 
frequency  of  the  boundary  disturbance  is  not  equal 
to  some  Fourier  mode  frequency),  one  can  expect  a 
smooth  flow,  desctibable  with  a  relatively  small  num¬ 
ber  of  grid  points  in  the  axial  direction  or  Fourier 
modes  if  a;=:0(l). 

The  results  in  the  present  work  are  for  nonresonant 


frequendes  u  =  1.0  and  u  =  2.5.  These  cases  are 
chosen  because  the  first  several  modes  dominate  the 
solution  and  only  the  first  six  modes  are  needed  to 
resolve  the  flow  field  in  the  former  case.  At  w  =  2.5, 
the  energy  is  more  equally  distributed  among  several 
modes,  sc  that  the  time  response  curves  (not  shown 
here)  are  more  irregular.  This  suggests  that  quite 
complex  time  response  curves  are  likely  to  be  found 
in  systems  containing  complex  acoustic  wave  fields. 

The  conceptual  approach  used  here  can  be  ex¬ 
tended  to  disturbances  driven  by  sidewall  injection 
transients,  reminiscent  of  burning  rate  variation  in 
real  systems.  Further,  it  is  likely  that  the  formula¬ 
tion  will  be  useful  for  a  reconsideration  of  flow  turning 
and  velodty  coupling  concepts.  Finally,  the  variable 
scaling  reported  here  will  enable  the  development  of 
more  accurate  numerical  methods. 
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Abstract 

A  mathematical  model  is  formulated  to  describe  the  initiation  and  evolution 
of  unsteady  vorticity  in  a  low  Mach  number  (M),  weakly  viscous  internal  flow 
sustsuned  by  mass  addition  through  the  side  wall  of  a  long,  nairrow  cylinder. 

An  0(M)  acoustical  disturbance,  generated  by  a  prescribed  harmonic  transient 
endwall  velocity,  interacts  with  the  basically  inviscid  rotational  steady  injected 
flow  to  generate  time  dependent  vorticity  at  the  side  wall.  The  steady  radial 
velocity  component  convects  the  vorticity  into  the  flow.  The  axial  velocity  asso- 
'  ciated  with  the  vorticity  field  varies  across  the  cylinder  radius  and  in  particular 
has  an  instantaneous  oscillatory  spatial  distribution  with  a  characteristic  wave 
length  0(M)  smaller  than  the  radius.  Weak  viscous  effects  cause  the  vorticity 
to  diffuse,  as  the  spatial  structures  are  convected  from  the  wall  toward  the  axis. 

The  magnitude  of  the  transient  vorticity  field  is  larger  by  0{M~^)  than  that 
in  the  steady  flow. 

An  initial-boimdary  value  formulation  is  employed  to  And  nonlinear  im- 
steady  solutions  when  a  pressure  node  exists  at  the  downstream  exit  of  the 
cylinder.  The  complete  velocity  consists  of  a  superposition  of  the  steady  flow, 
an  acoustic  (irrotational)  field  and  the  vortical  component,  all  of  the  same 
magnitude.  The  formulation  and  results  provide  a  conceptual  framwork  for  the 
study  of  injected  cylinder  flow  dynamics  which  supplements  traditional  acoustic 
stability  theory  by  providing  information  about  the  generation  and  evolution 
of  vortical  structures. 

‘Research  is  supported  by  the  Air  Force  OfRce  of  Scientific  Research  through  AFOSR  89*0023 
^Mailing  address:  Graduate  School  B-26;  University  of  Colorado  at  Boulder;  Boulder,  CO  80309 
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1  Introduction 


« 

* 


The  objective  of  this  study  is  to  demonstrate  how  and  where  unsteady  vorticity  is 
created  by  an  interaction  between  axially  propagating  acoustic  waves  and  a  flow  field 
caused  by  steady  sidewall  mass  injection  in  a  finite  length  cylinder.  Additionally, 
the  formulation  and  analysis  describes  the  transport  and  time-history  of  the  spatial 
distribution  of  vorticity  within  the  cylinder.  Finally,  the  results  show  that  transient 
rotational  flow  effects  are  crucial  to  the  evolution  and  stability  of  interned  fluid  dy¬ 
namics  when  the  characteristic  cylinder  Reynolds  number  (Re)  and  Mach  number 
(M)  are  very  large  emd  smadl,  respectively. 

The  mathematiced  model  describes  the  transient  response  of  a  flow  field  in  a 
cylinder  with  one  closed  end,  induced  by  strong  normal,  steady  sidewall  injection.  In 
particular  the  injected  flow  is  characterized  by  a  radial  speed  much  larger  than  the  ra¬ 
tio  of  the  characteristic  axial  speed  to  Re^^*.  Following  Taylor  (1956),  Culick  (1966a) 
shows  that  the  steady  flow  is  described  in  a  first  approximation  by  an  inviscid,  rota¬ 
tional  equation  system.  Acoustic  disturbances  of  0(M)  are  created  by  a  prescribed 
time-dependent  axial  velocity  variation  at  the  closed  end  of  the  cybnder.  They  prop¬ 
agate  through  the  basic  inviscid  shear  flow  field,  and  perhaps  unexpectedly,  create 
significant  vorticity  at  the  surface  of  the  porous  cylinder.  The  radial  component  of 
the  steady  injected  flow  field  carries  the  vorticity  into  the  entire  cylinder. 

Predictions  of  flow  dynamics  in  closely  related  systems  have  been  obtained  in 
the  past  from  traditional  acoustic  stability  theory  as  formulated  and  described  by 
Grad  (1949).  Culick  and  co-workers,  (See  Culick  (1990)  for  a  comprehensive  review) 
in  particular,  have  developed  since  the  mid-1970’s  an  increasingly  sophisticated  lin¬ 
ear  and  weakly  nonlinear  stability  theory  to  describe  disturbance  behavior  observed 
in  laboratory  e^eriments  and  practical  systems  like  solid  rocket  engine  chambers. 
Comparisons  with  experiments  appear  to  be  reasonably  good,  although  the  theory  as 
formulated  cannot  describe  the  generation  and  evolution  of  vorticity. 

Recent  experimental  (Brown  et  al.  1986a,b,  Brown  and  Shaeffer,  1992)  and  com¬ 
putational  (Vuillot  and  Avalon, 1991;  Vuillot,1991)  studies  include  definitive  evidence 
of  significant  transient  rotational  flow  structures  distributed  throughout  the  injected 
cylinder  flow.  In  addition.  Price  and  Flandro  (1993)  and  Flandro  and  Roach  (1992) 
have  formulated  approximate  mathematical  models  for  describing  the  generation  and 
evolution  of  vorticity  in  an  injected  cylinder  flow  field.  These  contributions  provide 
the  motivation  for  developing  a  comprehensive  flow  dynamics  model  which  includes 
both  the  irrotational  acoustical  phenomena  of  traditional  stability  theory  and  impor¬ 
tant  rotational  flow  effects. 

Stability  predictions  for  these  flow  systems  are  usually  obtained  from  mathemati¬ 
cal  models  that  reflect  linear  and  weakly  nonlinear ,  inviscid ,  irrotational  acoustic  con¬ 
cepts  (Ci^ck,  1966a, 1966b,  1967,1968,1970,1971,1973,1975,1976,1988,1990;  Williams, 
1985;  Price  and  Flandro,  1993;  Jahnke  and  Cuhck,  1993).  The  basic  acoustic  waves 
propagate  through  a  fluid  at  rest,  in  the  first  approximation,  and  are  described  in 
terms  of  Fourier  series  of  eigenfunctions  that  satisfy  totally  impermeable  boundary 
conditions.  Strictly  speaking  ,  these  solutions  do  not  accomodate  sidewall  injection 
nor  exit  plane  flow  at  the  downstream  end  of  the  cylinder.  As  a  result  the"  model 
cannot  account  for  an  interaction  between  the  acoustic  signals  and  a  sidewall  in¬ 
jected  cylinder  flow,  now  known  to  be  a  source  of  vorticity  generation  (  Flandro  and 
Roach,  1992). 

Flandro  (1974)  recognized  that  rotational  flow  effects  play  a  role  in  relatively  thin 
acoustic  boundary  layers  where  viscosity  is  of  significance.  Related  studies  for  inert 
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flows  have  been  carried  out  by  Tien  (1972),  Flandro  (1986),  Hegde  et  al.(1986)  and 
Price  and  Flandro  (1993).  Chemically  and  thermally  active  acoustic  boundary  layer 
flows  are  described  by  Hegde  and  Zinn  (1986),  Sankar  et  aJ.(  1988a, b),  Chen  et  al. 
(1990)  and  Matta  and  Zinn  (1993).  In  these  studies  the  boundary  layer  responds 
passively  to  externally  imposed  disturbances.  The  investigations  are  motivated  by 
a  need  to  understand  how  energy  is  exchanged  between  the  acoustic  disturbances 
and  mean  flow  as  fluid  injected  normally  from  the  wall  is  turned  towards  the  axial 
direction.  Until  recently,  conceptual  understanding  of  this  flow  turning  process  has 
been  based  largely  on  the  viscous  properties  of  the  thin  acoustic  boundary  layer. 

Signiflcant  efforts  have  been  made  to  develop  computational  models  for  acoustic 
boundary  layer  processes.  Baum  and  Levine(1987),  Baum  (1990),  Vuillot  and  Avalon 
(1991)  and  Vuillot  (1991)  have  employed  complete  Navier-Stokes  solvers  to  evaluate 
the  general  cylinder  flow  response  to  imposed  disturbances.  The  last  two  cited  works 
in  particular  demonstrate  that  rotationsd  disturbances  can  exist  in  a  substantial  por¬ 
tion  of  the  cylinder  for  appropriate  values  of  Re  and  M.  In  other  words,  the  vorticity 
distribution  is  not  always  confined  to  the  traditional,  viscous  acoustic  boundary  layer 
adjacent  to  the  cylinder  wsdl.  Rather,  for  appropriate  parameter  ranges  it  appears  to 
be  convected  from  the  wall,  out  into  the  main  portion  of  the  basically  inviscid  cylin¬ 
der  flow  by  the  mean  injected  flow  field.  Clearly,  the  pervasive  presence  of  vortical 
structures  in  the  intern^  flow  field  has  significant  consequences  for  the  conceptuad 
validity  of  traditional  irrotational  acoustic  stability  models. 

The  presence  of  rotational  inviscid  “acoustic”  disturbances  in  a  laboratory  flow 
field  was  first  implied  by  the  experiments  of  Brown  et  al.  ( 1986a, bl  and  Brown  and 
Shaeffer  (1992).  The  book  manuscript  of  Price  and  Flandro  (1993)  contains  an  im¬ 
portant  discussion  of  these  issues,  and  several  idealized  mathematical  models  are  used 
to  predict  the  structure  of  basically  inviscid  ,  but  weakly  viscous,  vortical  flow  struc¬ 
tures.  Mote  recently,  Flandro  and  Roach  (1992)  have  extended  the  theoretical  effort 
to  generate  a  purely  inviscid  model  for  boundary  generated  rotational  disturbances 
that  arise  from  an  inviscid  interaction  between  an  injected  flow  field  and  axial,  planar 
pressure  waves  propagating  in  the  cylinder.  Here  again  the  vortical  structures  are 
convected  out  into  the  cylinder  by  the  radial  component  of  the  injected  velocity  field. 
The  length  scsde  of  the  local  vorticity  variation  is  0(M)  smaller  than  the  cylinder 
radius,  implying  that  viscous  shear  stress  may  be  important  on  the  local  level.  This 
work  suggests  that  the  interaction  of  “acoustic  disturbances”  with  a  strongly  injected 
flow  field  leads  to  the  presence  of  important  vortical  structure  in  the  entire  flow  field. 

The  present  study  is  an  initial  step  in  formulating  a  rational  mathematical  model 
for  internal  flow  dynamics  which  incorporates  both  acoustic  phenomena  and  vortical 
flow  structures.  Perturbation  methods  are  used  to  derive  systematic  approximations 
to  the  complete  compressible  Navier-Stokes  equations.  An  initial- boundeiry  value 
approach  is  used  to  formulate  a  generalized  unsteady  mathematical  model  capable  of 
describing  both  non-resonant  and  resonant  time  history  of  solutions.  Fourier  series 
representations  of  the  velocity  and  pressure  are  obtained  in  terms  of  eigenfunctions 
that  satisfy  all  the  prescribed  boundary  conditions.  The  complete  aixial  velocity  is 
found  from  a  superposition  of  three  components  of  equal  magnituide.  First,  the  steady 
component  arises  from  a  solution  to  inviscid,  rotational  Euler  equations  known  by 
Culick  (1966a).  Secondly,  there  is  a  planar  irrotational  acoustic  field,  derived  from  a 
traditional  linear  wave  equation  which  satisfies  boundary  conditions  at  the  closed  and 
open  ends  of  the  cylinder.  Finally,  there  is  a  rotational,  weakly  viscous  component 
which  varies  on  two  distinctly  different  length  scales.  The  global  spatial  distribution 
occurs  on  the  radial  length  scale,  while  loc^y  there  are  velocity  gradients  on  a  scale 
0(M)  smaller.  As  a  result  the  vorticity  of  the  transient  field  is  0{M~^)  larger  than 
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that  in  the  steady  field.  The  convection  of  the  orticity  is  found  to  be  describea  by 
a  linear  equation.  Most  significantly,  the  diffusion  of  vorticity  on  the  shorter  length 
scale  is  described  by  a  nonlinear  di^usion  equation.  Solutions  are  found  numerically 
by  using  a  finite  difference  method,  as  well  as  semi- analytically  with  a  Galerkin-like 
method  developed  by  Wang  and  Kassoy  (1990a, b,c,  1992a, b,  1993). 

Evaduation  and  interpretation  of  the  full  results  show  that  rather  complex  vortical 
structures  axe  present  throughtout  the  cylinder  sufficiently  long  after  the  disturbance 
is  initiated  at  the  end  wall.  When  the  input  from  several  different  Fourier  modes  is 
important,  the  time  response  curves  are  quite  irregular  and  may  be  mistaken  for  a 
“turbulent”  response. 

Fully  computationad  methods  are  used  by  Kirkkopru  et  al.(1994)  to  provide  sup¬ 
porting  evidence  for  the  solutions  found  here  by  quasianalytical  means. 


2  Mathematical  Formulation 


A  cylindrical  tube  with  side  wall  injection  is  shown  schematically  in  Figure  1.  in 
which  L  is  the  axial  length  and  D'  is  the  diameter.  The  radius  is  /?’  =  D'  12.  The 
complete  non-dimensional  equations  describing  the  fluid  dynamics  and  acoustics  for 
an  axisymmetric  system  can  be  written  in  the  form 
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and  $  is  the  viscous  dissipation  function.  The  non-dimensional  variables  are  defined 
in  terms  of  dimensional  quantities  (with  a  prime)  by 
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The  reference  vedue  pj,  measures  the  bcisic  static  pressure  m  the  cylinder,  while  the 
analogous  density  and  temperature  valuer  p^,  Tg  respectively  represent  properties 
of  the  injected  fluid.  The  known  characteristic  injection  speed  V'^g  is  related  to  the 
derived  characteristic  axial  speed  V’g  by  the  approximate  mass  conservation  relation¬ 
ship  where  the  aspect  ratio  ^  Characteristic  length  scales  for  the 

axi2d  and  radial  dimensions  are  defined  by  L'  and  R'  respectively.  Time  is  nondimen- 
sionalized  with  respect  to  the  axial  acoustic  time  t'  =  where  Cg  =  {'^R'Tq)^  is 

the  characteristic  sound  speed.  The  reference  material  properties  fcg,  pg  and  C„o  are 
defined  at  temperature  Tg.  The  parameter  7  is  the  ratio  of  specific  heat  and 

D-  _  p  _  Po^po  _Yi£.  t’r\ 

^0  ^0 
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where  typically  the  Prandtl  number  Pr=  0(1),  the  axial  Mach  number  M  <g:  1  and 
axiad  Reynolds  number  Re  1. 

Initi^ly,  a  steady  flow  exists  in  the  system,  driven  by  spatially  distributed  normal 
injection  from  the  wall  where  the  no-shp  condition  is  satisfied.  Symmetry  prevails 
along  the  axis  at  the  cylinder.  The  closed  end  of  the  cylinder  is  impermeable  and  for 
simplicity  the  exit  is  assumed  to  be  a  pressure  node.  The  mathematical  form  of  the 
boundary  conditions  may  be  written  as 

dV 

r  =  0;  V;  =  =  0,  r  =  1;  V;  =  -KM,  K  =  0,  T  ^  1  (8) 

2  =  0;  V,  =  0,  2  =  1;  P  =  l.  (9) 

The  steady  flow  is  disturbed  at  2  =  0  by  imposing  a  harmonic  end  wall  axial 
velocity  variation  that  is  independent  of  the  radial  coordinate; 

2  =  0;  V’,  =  i4sinu>t,  t  >  0;  0  <  r  <  1,  (10) 

where  the  amplitude  A=0(1). 

It  should  be  noted  that  the  imposed  end  wall  disturbance,  of  the  same  order 
of  magnitude  as  that  of  the  steady  axial  speed,  induces  both  mechanical  and  ther¬ 
modynamical  disturbances  of  similar  magnitude  into  the  gas.  These  relatively  large 
variations  cannot  be  described  by  linear  theory  alone.  Therefore,  the  weakly  nonlinear 
theory  developed  ^.ere  differs  from  traditional  small  disturbance  theory. 

Strong  injection  prevails,  in  the  sense  that  V^g  »  (Cole  and  Aroesty  1967), 

Rei 

which  implies  that  Re  P.  Also  note  that  in  this  modelling  effort  it  is  assumed 
that  5  1,  to  describe  a  long,  nairrow  cylinder. 


3  Steady  State  Flow 

The  steady  state  flow  generated  by  time  independent  mass  addition  on  the  side  wall 
can  be  described  in  terms  of  the  variables: 


{P,p,T)  ~  1  -f  M^{Po$, poi,Tot)  + 


(v;,v;)~(v;o.,v;o.)  +  ---;  (n) 
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valid  for  the  asymptotic  limit  A/  — »  0.  Equation  (11)  can  be  used  in  (1) — (5)  to  find 
the  leading  order  equations  describing  an  incompressible,  inviscid,  rotational  flow  that 
satisfies  the  no-slip  and  injection  boundary  conditions  on  the  side  wall  and  symmetry 
conditions  on  the  axis,  given  in  (8)  and  (9) 


1  d{rKo,)  dKo,  ^ 

r  dr  dz 

Po,  =  Po,{z) 

,,  dV,o,  ,  „  dV,o.  1  dPo. 
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(12) 

(13) 

(14) 


The  transport  terms  are  excluded  from  the  leading  order  equations  because  ^  1. 

Equation  (13)  arises  because  the  aspect  ratio  5  »  1.  The  solutions  for  the  radial  and 
axial  velocity,  as  well  as  the  pressure  distribution  can  be  written  in  the  form 


Vm,  =  (irjT  V;w(’’)i(t)  cos(^r’), 


(15) 

(16) 


where  —Vr„{z)  is  an  arbitary  time-independent  side  wall  injection  distribution.  Re¬ 
lated  solutions  can  be  found  in  Culick  (1966a),  Price  and  Flandro  (1993)  and  Flandro 
and  Roach(1992).  It  should  be  noted  that  recently  Balakrishnan  et  al.(1991)  derived 
a  related  solution  valid  for  Af  =  0(1)  <  1. 


4  Core/Boundary  Layer  Solutions 

For  certain  parameter  values  the  flow  response  to  the  endwall  driven  disturbance  in 
(10)  can  be  described  in  terms  of  an  irrotational  acoustical  core  and  a  relative  thinner 
weakly  viscous,  rotational  boundeiry  layer  adjacent  to  the  cylinder  wall.  Aspects  of 
this  traditional  model  have  been  described  by  Wang  and  Kassoy(  1992a).  Here,  the 
objective  is  to  use  the  solutions  themselves  to  determine  the  conditions  for  which  the 
core-boundary  layer  structure  is  valid. 

The  asymptotic  expansions  for  the  unsteady  core  flow  can  be  written  as 

+  (v;,v;)~EJV"(v„,v'„)i  (i7) 

n=0  n=0 

in  the  limit  M  —*  0.  Equation  (17)  can  be  used  in  (1) — (5)  to  derive  the  lowest  order 
acoustic  system  valid  in  the  limit  M  —*  0,  Re  —*  oo,  6M  =  0(1); 
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The  velocity  components  consist  in  general  of  both  a  steady  state  part  and  an 

unsteady  part  respectively;  (Ko,Ko)  =  (KojiKo*)  +  (Ko>Ko)  •  It  follows  that  the 
unsteady  flow  equations  can  be  derived  by  substracting  the  steady  state  equations 
(12)-(14)  from  (18) — (21).  Given  the  boundary  condition  in  (10)  the  lowest  order 

radial  speed  Vro  =  0. 


4.1  Linear  Planar  Solution 

The  unsteady  part  of  the  leading  order  equations  can  be  combined  into  a  planar  wave 
equation  in  terms  of  the  axial  velocity  component; 


at*  dz^ 

The  corresponding  initial  and  boundary  conditions  are  given  by 


(22) 


II 

o 

o 

II 

o 

o 

II 

o 

(23) 

2  =  0, 

Vgo  =  sinu;t; 

(24) 

2  =  1, 

(25) 

The  simplicity  of  the  equation  can  be  attributed  to  the  large  aspect  ratio  condition 
^  3>  1  .  The  boundary  condition  at  the  exit  is  derived  from  the  pressure  node 
condition  in  (9). 

The  general  solution  for  V^o  is  ; 

Vto{2,t)  =  sinwt 


-  I  sin(6„.t)  +  tcos(h„.t)|  sin(6„.z), 


(26) 


where  =  (n  +  |)ir,  and  the  last  set  of  terms  represents  a  resonant  effect  present 
only  when  u;  =  6„..  A  careful  look  into  the  solution  provides  us  with  some  insight 
into  the  properties  of  the  acoustical  flow. 

•  The  first  term  itself  and  the  second  part  of  the  nonresonant  Fourier  series  repre¬ 
sent  quasi-steady  motion  oscillating  at  the  driving  frequency.  The  other  Fourier 
series  terms  can  be  decomposed  into  two  counter-propagating  planar  travelling 
waves. 


•  If  the  driving  frequency  a;  is  not  equal  to  one  of  the  natural  frequencies  6„,  the 
solution  is  bounded.  If  u;  is  very  close  to  one  of  these  natural  frequencies,  then 
beats  will  appear  due  to  the  interaction  between  the  quasi-steady  motion  and 
one  pair  of  travelling  waves. 
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•  Resonance  occurs  when  a;  =  6„,  and  the  amplitude  of  one  mode  grows  linearly 
with  time. 

Table  1  contains  results  for  a  system  where  =  10~^  s  so  that  dimensional 
frequencies  can  be  considered.  When  u;  159  Hz,  the  response  shown  in  Figure  2 
for  Ko  at  2  =  0.5  is  bounded  and  the  contributions  are  mainly  from  the  first  harmonic 
oscillation  term,  the  first  few  axial  modes  and  quasi-steady  modes  in  (26).  Beats  are 
observed  in  Figure  3  when  uj  ^  238  Hz.  The  period  of  the  beat  is  about  90s.  due 
to  the  interaction  between  the  quasi-steady  modes  and  the  first  axial  mode.  Linear 
monotonic  amplitude  growth  seen  in  Figure  4  is  primarily  from  the  resonant  axial 
mode  in  (26)  when  u  =  250  Hz. 

The  pressure  solution  Pq{z,  t)  can  be  obtained  from  a  first  integral  of  the  unsteady 
part  of  (18)  and  the  isentropic  relationship  Pq  =  7po- 

4.2  Boundary  Layer  Correction 

The  leading  order  core  acoustic  solution  does  not  satisfy  the  no- slip  boundary  con¬ 
dition.  Under  certain  conditions  to  be  defined  quantitatively,  the  transition  to  zero 
axial  velocity  at  the  wall  occurs  in  a  relatively  thin  boundary  layer  which  has  a  multi¬ 
ple  scale  structure.  In  particular  the  overall  radial  thickness  of  the  layer  is  defined  by 
viscous  considerations.  But  within  it  there  is  a  smaller  length  scale  associated  with 

the  distance  traveled  by  an  injected  fluid  particle  on  the  time  scale  Due 

to  the  relatively  large  injection  condition  »  U^'o/iZea  the  boundary  layer  flow 
is  inviscid  and  rotational  in  the  first  approximation.  Viscous  stresses  appear  in  the 
second  order  description,  but  are  essential  to  finding  the  complete  solution,  as  might 
be  expected  in  a  multiple  scale  analysis. 

The  multiple  scale  structure  is  defined  in  terms  of  the  variables; 

1  —  r  1  —  r 

where  0  =  ReM^fS^  and  M  1  when  the  core/boundary  layer  concept  is  valid. 

The  partial  derivatives  with  respect  to  r  must  be  replaced  by 

(i)-i(a- 

(I?)  (af^)  + 

The  variables  are  expanded  as 

(,P,p,T)  '  1  +  E  (30) 

n3=0 

V,'-Ko  +  jV,i+---,  v;  -  +  • .  • ,  (31) 


8 


The  lowest  order  boundary  layer  equations  are; 

dV,o  ,  ,,  dV^  1  dPo 

dt  ^  5^  -i  dz' 

Po  =  Po{z,t). 


(32) 

(33) 


Equation  (32)  is  an  inviscid  rotational  equation  which  can  satisfy  the  no-slip  boundary 
conditions  oa  ^  =  t}  =  0.  Vorticity  transport  can  be  described  by  the  first  derivative 
of  (33)  with  respect  to  keeping  in  mind  that  and  ^  are  independent  of  the 
short  length  scale  variable.  The  resulting  equation  shows  that  vorticity  is  convected 
invariantly  by  the  radial  wall  injection  velocity  into  the  boundary  layer  along  well 
defined  characteristic  lines,  (p  =  t  —  (^)- 

The  second  order  momentum  equation  is  obtained  from  terms  of  0(y ); 


dv^i  ^  ^  dK,  _  dv^o  ^  a^v.0 
dt  dri  de 


(34) 


where  a  viscous  stress  term  is  present  and  the  pressure  gradient  term  is  absent  since 
^  >  A/. 

The  acoustic  solution  in  the  core  provides  an  outer  boundary  condition  and  the  no¬ 
slip  condition  on  the  side  wall  provides  an  inner  boundary  condition  for  (32)  smd  (34). 
The  leading  order  acoustic  core  solution  shows  that  all  the  terms  can  be  classified  into 
the  following  two  forms:  Ko(«)c*”‘  with  fl  =  u;  or  6„,  and  -t  cos(6„.<)  sin(6„.2).  Thus 
the  boundary  conditions  are 


(  =  V  =  0,  Ko  =  0;  (35) 

^,7} -*  00-,  V;o~Ko(2,0-  (36) 

Equations  (^32) — (34)  and  (35) — (36)  are  used  to  find  quasi-steady  solutions  to 
every  driving  irequency  in  the  core  solution.  For  any  of  the  nonresonant  modes 
n  =  u>  or  &„  but  uj  ^  bn  for  any  integer  n,  the  boundary  layer  solutions  can  be  written 
as 


V,o  =  f  (f ,  <1.  V.,  =  G({, ,,  (37) 

These  solution  forms  can  be  substituted  into  (32)  and  (34)  to  determine  F  and  G. 
The  former  found  from  (32)  ,  (35)  and  (36),  is 

=  +  (38) 

Vrw 

where  the  undetermined  coefficient  function  C{ti,  z)  must  satisfy  the  conditions 


=  0;  C  =  -Ko(«). 

Tj  — ►  oo;  (7  =  0. 

Equation  (34)  can  then  be  rewritten  in  terms  of  G  ajid  C  as 

in 


5G  ^  in  _  . 


dC  JP 


(39) 

(40) 


(41) 
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In  order  to  avoid  secular  growth  of  G  with  respect  to  the  variable  the  quantities  in 
the  square  bracket  must  be  set  to  zero.  Therefore, 


, 

together  with  boundary  conditions  (39)-(40)  are  solved  to  find 


C{v,z)  =  -y,oexp(- 


VI. 


v)- 


The  multiple  scale  solution  for  a  given  single  frequency  fi  has  the  form: 


V;o(^,  Vy  0  =  -Vzo{z)  \  exp 


n*  in  ^ 
I  KIW’' 


1-1 


.iOt 


(42) 


(43) 


(44) 


where  — V^ui(^)  is  the  steady  side  wall  injection  velocity.  The  exponential  term  rep¬ 
resents  the  vortical  part  of  the  axial  unsteady  flow  in  the  boundary  layer.  The  com¬ 
bination  of  this  vortical  component  and  e*”*  term  can  be  written  as: 


exp 


in 


^  -)-  int 


(45) 


where  <fi  =  t  -  ^  is  the  constant  phase  line  or  characteristic  line  for  the  unsteady 
axial  vortical  velocity.  The  radial  traveling  speed  for  constant  ^  line  can  be  therefore 
described  by: 


dr 

Tt 


(46) 


This  shows  explicitly  that  the  vorticity  is  convected  by  the  steady  injected  fluid 
velocity  in  the  boundary  layer. 

The  first  part  of  the  exponential  term  in  (45)  describes  amplitude  damping  arising 
from  viscous  effects  because  the  17- variable  denned  in  (27)  is  scaled  with  respect  to 
Re,  in  part.  The  second  part  describes  harmonic  spatisd  oscillations. 

When  resonant  driving  is  present  n  =u/  =  bn-,  and  has  the  form 


VzoiLny^yt)  =  -tsin(6n-»)sin(h„.t) 

+  {t  sin(fc4)sin(6„*t)  +  5n*t  cos(fc^)  cos(6„t) 


^  cos(k()  +  ^ 


I 

u 


(1  -I-  bn-) 

Vriiz) 


sin(fc^) 


sin(5n*0' 


2(1  -t-  bn-) 


2k^{l  +  bn-) 

(  cos{k()  -  -  7:^^)  sin(fcO 


1  6* 
cos{bn-t)}  —  sin(6n*2)exp(— r;r2^77) 


(47 
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wliere  k  = 

When  (  =  rj  =  0,  the  solutions  satisfy  the  no-slip  boundary  condition  on  the  wall. 
On  the  other  hand,  when  ^  and  ?/  — »  oo,  the  core  solution  is  recovered  in  an  oscillatory 
manner  since  the  amplitude  of  the  exponential  term  goes  to  zero  harmonically.  The 
effective  thickness  of  the  boundary  layer  depends  strongly  on  Q  and  V^.  A  large  value 
of  n  promotes  relatively  rapid  exponentii  decay,  implying  that  a  high  frequency 
disturbance  is  associated  with  a  thinner  transition  boundary  layer.  Alternatively,  low 
frequency  forcing  fosters  thick  boundary  layers.  Thus,  higher  order  modes  tend  to 
be  associated  with  effectively  thinner  boundary  layers.  The  same  type  of  argument 
demonstrates  that  increasing  the  value  of  Vu,{z)  enhances  the  overall  boundary  layer 
thickness. 

A  complete  solution  for  the  axial  velocity  in  the  boundary  layer  consists  of  an 
infinite  sum  of  terms  obtained  from  (44)  amd  (47)  ;  one  for  each  frequency  u;  and 
6„  in  (26).  The  spatial  structure  of  such  a  solution  will  be  quite  complex,  given  the 
oscillatory  dependence  on  the  value  of  It  is  perhaps  more  illustrative  to  look  at 
the  results  for  a  single  frequency. 

The  reduced  axial  velocity  inside  the  vortical  layer,  is  plotted  against  ( 

in  Figure  5  with  A/  =  0.01,/3  =  0.1  and  =  1  for  f2  =  2.5  and  Q  =  3.0.  The  core 
solution  is  recovered  at  about  ^  =  10  for  fl  =  2.5  which  correspondes  to  r  =  0.9  . 
In  contrast,  the  boundary  layer  thickness  is  a  little  smaller  for  the  higher  frequency 
n  =  3.0.  Of  course  the  overall  boundary  thickness  is  determined  by  the  lowest  mode 
in  the  system.  The  core/boundary  layer  structure  is  valid  only  when  is 

sufficiently  small,  where  <  1,  to  permit  exponential  decay  in  (45)  to  occur  close  to 
the  cylinder  wall  (say  r  >  0.9).  1(0  =  0(1),  and  is  sufficiently  large  and/or  Cl  is 
sufficiently  small,  then  a  new  multiple  scale  perturbation  technique  is  needed  to  find 
solutions  where  rotational  effects  fill  the  entire  chamber. 


5  Rotational  Flow  Across  the  Cylinder 

5.1  Problem  statement 

The  objective  now  is  to  develop  a  mathematical  model  for  the  disturbance  evolution 
when  rotationsJ  effects  are  significant  across  the  entire  cylinder.  The  wymptotic 
expansions  for  the  velocity  components  and  thermodynamic  variables  in  the  limit 
Af  — »  0  are 


V,  ~  V,o,{z,r)  +  Ar'V^{z,r,t)  (48) 

n=sO 

Vr  -  Vro,{z,r)  +  M^Vrn{z,r,t)  (49) 

n=!0 

{P,P,T)  ~  l  +  M'£M^{Pn,Pn,9n)  (50) 

n=0 

It  is  recognized  that  two  disparate  length  scades  are  important;  the  tube  radius 
and  a  much  shorter  length  associated  witn  the  radial  distance  traveled  by  a  fluid 
particle  on  the  acoustic  timescale.  A  multiple  scale  analysis  will  be  carried  out  in 
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terms  of  the  variables  ri  and  r2  defined  by 


ri  =  l-r;  rj  =  /  — t-  da.  (51) 

Jo  -M\ro,{<r) 

The  second  transformation  includes  an  integral  of  the  steady  radial  velocity  field  for 
the  case  of  constant  steady  wall  injection  V;„  =  1  which  simplifies  the  describing 
equations  considerably.  It  is  noted  that  when  the  center  line  is  approached  rj  — ►  1 
the  integral  diverges  and  — ►  oo. 

Each  of  the  dependent  variables  is  written  in  terms  of  ri  and  r2  instead  of  r  alone. 
The  partial  derivatives  with  respect  to  r  in  equations  (l)-(5)  must  be  replaced  by 


5.2  Lowest  order  mathematical  model 

The  relations  (48) — (50)  can  be  substituted  into  (1) — (5)  to  find  the  leading  order 
equations  in  the  umit  M  —*  Q.  First,  the  spatially  nomogeneous  boundary  forcing  in 
(10)  and  the  condition  5  >  1  imply  that  Ko  =  0.  Then, 


dRo 

.dRo 

dV^  f  1  \  dVri 

(54) 

(55) 

dt 

dr. 

dz  WroJ  dr,  ’ 

dV,o  dVa  IdPo 

dt  dr,  7  dz  ’ 

dPo  dPo 
dri  ~  dr,  ~  ' 

(56) 

d9o  dOo  _  (7  -  1)  dPo 
dt  dr,  7  dt  ’ 

(57) 

Pq  =  Ro  +  9o- 

(58) 

Following  a  procedure  related  to  that  described  by  Lagerstrom(1964),  and  similar 
to  that  employed  by  Price  and  Flandro  (1993)  and  Flandro  and  Roach  (1992),  the 
variables,  except  for  Pqi  &re  divided  into  irrotational  planar  and  rotational  nonplanar 
parts, 

Vxo  -  Wop{z,  t)  +  Wo{z,  t,  ri,  Vi), 

Ro  =  Rqp{z,  t)  +  Ro{z,  t,  ri,  r2), 

^0  =  ^op(«,0  +  ^o(«,t,ri,r2). 
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(59) 

(60) 
(61) 


Equations  (59) — (61)  can  be  used  in  (54) — (58)  to  show  that  the  planar  functions  are 
described  by 

-  —IT-  <«■*> 

Po  =  Rop  +  ^op(  (65) 

a  system  nearly  identical  to  (18) — (21).  The  initial/boundary  conditions  are: 

t  =  0,  Wop  =  0,  -^=0,  (66) 

2  =  0,  lVop  =  sinu>t;  (67) 

^  =  1,  ^=0.  (68) 


2  =  1, 


in  analogy  to  (23) — (25). 

The  nonresonant  acoustic  solution  for  the  planar  contribution  is 

Wop{t,  z)  =  --52  a„  sin(4i;t)  -  sin(A„t))  8in(An2),  A„  ^  u;,  (69) 

7  n=l  V  / 

00 

^o(t,2)  =  2  an  (-cos(a;t)  +  cos(AnO)co8(A„z),  An  #  u;,  (70) 

nssl 

Po  =  IfRop  (71) 

where  Un  =  —  An  =  (n  —  ^)7r.  Equation  (69)  is  equivedent  to  (26)  for  the 

nonresonant  case.  The  first  term  in  the  sums  of  (69)  and  (70)  arises  from  forcing 
at  frequency  u,  and  the  second  term  describes  the  eigenfunction  response.  Only  the 
nonresonant  case  will  be  considered  in  the  present  work. 

The  equations  for  the  rotational  components  are 

(70\ 

dt  dr,  dz  \VroJ  dr,  ’  ^  ’ 

^  +  9o  —  (^5) 

Equations  (74)  and  (75)  can  be  combined  to  show  that  the  leading  order  vorticad  flow 

is  incompressible: 


,  dRo  _  ^ 

dt  ^  dr,  ~ 
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Therefore,  (72)  can  be  rewritten  as 


W  /  1  \ 

dz  VKoJ  dT2  ’ 

(77) 

which  can  be  used  to  find  Kj  once  Wq  is  known. 

The  relevant  initial  and  boundary  conditions  are 

t  =  0;  lVo  =  0,  ^=0, 

(78) 

N 

II 

_o 

II 

o 

(79) 

1  n 

^  =  1;  -^=0. 

(80) 

,  dWo  .  „ 

r,  =  l,  =  0,  ^^=0. 

(81) 

n  =  rj  =  0;  Wq  =  -  Wop{t,  z),  Sq  = 

(82) 

The  first  of  (81)  can  be  combined  with  (73)  and  the  initial  condition  (78)  to  prove 
that  =  0  on  the  axis  ri  =  0  for  all  t.  Equation  (82)  corresponds  to  the  no-slip 
condition  and  isothermal  flow  injection.  Equations  (73),  (74)  and  (76)  show  that  Wq, 
and  Rq  are  invariant  on  a  characteristic  line  defined  by 

r}  =  t-r2,  (83) 

but  vary  across  the  ij  lines.  The  tj  =  0  line  enters  the  cylinder  at  i  =  0  through  side 
wall  (rj  -  0)  and  subsequently,  aX  t  =  c,  t)  =  c  appears  at  r2  =  0.  At  a  particular 
time  T,  constant  t;  lines,  which  range  from  0  to  T,  are  transported  toward  the  zods 
by  convection  at  the  local  radial  steady  velocity,  as  found  from  a  time  derivative  of 
(83)  sdter  using  (51). 

The  invisdd  equation  in  (73)  can  be  combined  with  the  first  of  (82)  to  show  that 
the  wall  is  a  source  of  vorticity  because 

5W’o(t,x,0,0)  dWo  _dWop{t,z) 

dr2  ~  dt  dt  ,  ^  ^ 

where  Wop  is  known  from  (69).  It  is  noted  that  the  largest  unsteady  non  dimensional 
vorticity  term  is  given  by  f2s  =  The  parameter,  arises  from  large 

gradients  occuring  in  the  spatially  oscillatory  velocity  profile  on  the  short  length  scale 
r2.  Equation  (73)  also  shows  that  the  vorticity  generated  at  the  wall  is  convected  into 
the  cylinder  by  the  steady  radial  velocity  field  Vro,{r). 
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In  general,  solutions  to  the  inviscid  first  order  hyperbolic  equations  in  (73)-(75) 
can  be  written  symbobcally  as, 


1^0  =  n,  2),  V  =  i-  ^2  (85) 

^0  =  (86) 

Explicit  functional  dependence  of  the  variables  can  be  found  from  consideration  of  the 

higher  order  equations,  as  is  typicail  in  a  multivariable  study.  Once  Wq  is  found,  then 
the  mass  conservation  equation  (77)  can  be  integrated  with  respect  to  r2  to  find  the 
radial  velocity  V',i.  The  vortical  temperature  and  density  fields  can  be  found  using 
related  methods. 


5.3  Higher  order  consideration 

Equations  (48) — (50)  can  be  combined  with  (1) — (5)  to  find  the  0(M)  equation  set 
in  the  limit  M  0.  The  same  procedure  used  to  find  the  leading  order  solution  is 
employed  so  that  the  variables,  except  for  Pi,  are  divided  into  irrotational  planar  ^md 
rotational  nonplanar  parts, 


Vzi  =  V^ip{2,  t)  +  Wi{z,  t,  ri,  rj), 

Ri  =  Ripiz,  t)  +  Pi(z,  t,  ri,  ra), 
^1  =  Bip{z,t)  +  ^i(z,t,T*i,r2). 

The  planar,  acoustic  equations, 


dR\p 

dt 

dt 

dt 

Pi 


^  +  R^Wop] , 


^  f  l,p  p  1  ^ 

dt 


(7  -  i)^ip  + 


(7  -  2)(7  -  1)  p2 

0 


Rip  +  ^ip  +  Ritp^op 


(87) 

(88) 
(89) 


(90) 

(91) 

(92) 

(93) 


containing  quadratic  driving  terms  associated  with  lower  order  acoustics  are  not  con¬ 
sidered  furtW  here,  although  they  may  be  important  for  studying  acoustic  streaming 
effects. 

The  largest  possible  viscous  effect  occurs  when  5*/Re=M*,  in  which  case  the 
higher  order  axial  momentum  equation  for  Wi  has  the  form 


dWi  dWi 

dt  ^  dT2 


VX  drl 


d^Wo,li,dPo  ,  ,,  ,d{VX  +  \\o)  ■ 

- ~  ^  ^ - ^z - 


+v; 


<0( 


dv^o,  ^ ,,,  dWop  ^ divx  +  v;o) 

+  ^Op  Q-  +  VrO»- 


dz 

Ki  \  dV,o 
dr2 


dz 


dr  I 


\VrOsJ 


-K 


av. 


xOi 


rOt' 


dvi 


(94) 
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the  analogue  to  (73).  It  is  (94)  that  will  give  us  information  about  the  behaviour  of 
the  leading  order  axial  speed  solution. 

The  corresponding  energy  equation  for  analogous  to  (74)  contains  a  conduction 
term.  Thus,  transport  effects  are  important  conceptually  in  the  distribution  and 
evolution  of  the  rotational  variables. 

If  the  transformation  of  the  coordinate  system  from  (t.ri.rj.z)  to  {rj,ri,r2,z)  is 
made,  then  the  derivatives  with  respect  to  t  and  t2  must  oe  replaced  by 


It  follows  that  (94)  can  be  written  as 


dWi 

dr2 


n 


1  d^Wo  a(Ko.  +  v^o) 

v*0 - 


VX  dri^ 

”  dz 


K 


rOt 


dz 

^(KlWo)  d{Vr,Wo) 


dr2 


dr} 


y  -  dWo^  dWo  I  ^  dPo 

02  Oz  Oti  7  Oz 


(96) 


An  integration  of  (96)  with  respect  to  r2,  holding  r},  rj  and  z  fixed  will  generate 
secular  growth  in  rj  umess  certain  terms  are  suppress^.  In  considering  the  impact  of 
each  term,  it  is  important  to  remember  that  the  harmonic  t-dependence  of  the  planar 
acoustic  solutions  in  (69)-(71)  must  be  rewritten  in  terms  of  t}  and  r2  by  using  the 
transformation  7  =  t  —  r2  in  (83).  When  written  in  the  coordinate  system  (z,  t,  rj,  r2) 
the  suppressed  terms  take  the  form; 


1  a-Wo  ,  „  awo  „,i,  BWo  ,i,  Bv,o.  „  aw^  „ 

+  ^rOf— - 2VVo— 5 - Wo— - V,0»-5—  =  0. 


VX  drl 


dri 


dz 


dz 


dz 


(97) 


which  is  a  nonlinear  diffusion  equation  for  the  rotational  axial  velocity  Wq  with  a 
time- like  variable  ri.  The  solution  for  Wo  must  satisfy  an  “initial”  condition  from 
(82) 


Wo(t,2,ri  =  O.rj)  =  -Wop{r},z)  77  >  0 

1  “  /A  \ 

=  -  £  I  —  sin(u;T7)  -  siniXnV)  sin(A„z),  A„  ^  u, 

7  \  ‘^  / 

Wq(t,2,ri  =  0,r2)  =  0  77  <  0  (98) 

and  a  boundary  condition  at  the  center  line  from  (81) 

n  — »,  ^  =  0.  '  (99) 


In  addition  a  condition  must  be  specified  on  ri  >  0,r2  =  0  which  is  compatable  with 
(98)  at  the  point  ri  =  r2  =  0.  This  is  necessary  because  ri  and  r2  are  treated  a« 
independent  variables  in  (102).  The  reasonable  choice  is  given  by 


Wo{t,Z,T2,r2=0)  =  -Wop{t,z). 


(100) 
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The  nonlinear  term  in  (97),  is*  i^resent  in  our  problem  because  the  0{M) 

boundary  disturbance  is  larger  than  that  used  in  earlier,  basically  linear  studies  (Flan- 
dro  and  Roach,  1993).  Its  presence  suggests  that  wave  steepening  and  other  forms  of 
instability  may  occur  in  the  evolving  flow  field.  If  the  imposed  endwall  disturbance  is 
smaller, and/or  axial  variations  are  ignored,  then  a  linear,  viscous  diffusion  equation 
is  derived,  which  is  related  to  that  used  by  Price  and  Flandro  (1992). 

One  solution  approach  is  based  on  a  finite  difference  approximation  to  (97)-(100). 

Solutions  for  Wo{z,  t,  rj,  rj)  can  be  found  by  employing  a  second  order  accurate  Adam- 
Bashforth/Crank-Nicolson  scheme.  Most  of  the  results  presented  here  (see  section  6) 
have  been  found  in  this  way. 

Given  the  forcing  condition  in  (98),  it  is  also  possible  to  find  solutions  in  terms  of 
the  eigenfunction  set  {sin(A„z)},n  =  1,2,  •  ■ 

OO 

^oit,z,rur2)  =  ^  An(t,ri,r2)sin(A„r)  (101) 

n=l 

Coupled  partial  differential  equations  for  the  Fourier  coefficients  A„,  are  found  by 
using  (101)  in  (97)  and  invoking  orthogonality  conditions  for  the  eigenfunction  set  on 
the  intervsd  [0,1].  The  results  are; 


2  /  ~ 

9ri  VrO,  dz  "  Vro, 

®nni  nj  ^  ^ni  ^  ~  ^ 

and  are  subject  to  conditions  obtained  from  (98)  and  (99). 

(a)  Initial  Condition: 

A.n{t,  ri  =  0,  rj)  =  ^a„  sina;(f  -  r2)  -  sin  An(f  -  r2) j 


=  0  r2  >  f 

(b)  Boundary  Conditions: 


r,  — ►  OO,  =  0 

ar2 


where 


>ln(i,ri,r2  =  0)  =  ia„  sin(u/t)  -  sin(A„t)j 

Onnin,  =2A„,  [  sin(A„, z)  cos(A„, z)  sin( A„r)d2 
Jo 

bnni  =  f  zcos(An,z)sin(Anz)dz 
Jo 


(102) 

0  <  r2  <  < 

(103) 

(104) 

(105) 

(106) 
'  (107) 


Equation  (102)  is  a  coupled  system  of  quasi-linear  diffusion  equations  with  a  time  like 
variable  r-i,  and  nonlinear  source  terms.  The  physical  time  t  is  a  parameter  of  the 
differential  equations,  appearing  explicitly  only  in  the  boundary  conditions.  Solutions 
to  (102)-(107)  are  described  in  Section  7-  A  system  truncation  approach  is  used  to 
find  a  finite  number  of  An’s. 
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6  Finite  Difference  Solutions  for  the  Nonlinear 
System 

Solutions  to  (97)-(100)  for  Wo{z,t,ri,ri)  have  been  found  by  using  a  finite  differ¬ 
ence  method  based  on  the  second  order  accurate  0(^rj,  5r|)  Adam-Bashforth/Crank 
Nicolson  scheme.  This  scheme  is  a  semi-implicit  amd  neutrally  unstable  in  the  sense 
that  it  works  well  if  the  nonlinear  effect  is  moderate  in  comparison  with  the  viscous 
effect.  For  example,  instability  results  if  the  amplitude  of  the  initial  condition  for  a 
given  ijj  is  larger  than  a  threshold  value.  The  real  solution  lies  on  the  locus  defined 
by  (51)  relating  rj  and  rj.  In  general  this  curve  does  not  coincide  with  the  compu¬ 
tational  grid  points.  Hence  cubic  spline  interpolations  in  r2  are  employed  at  every 
chosen  rj  value  to  obtain  the  desired  real  solution.  The  far  end  boundary  condition 
for  r2  — »  oo  is  implemented  by  providing  a  sufficiently  large  number  of  grid  points 
between  the  convected  outer  edge  of  the  rotational  layer  ajid  the  finite  location  of 
the  computational  boundary  with  respect  to  r2.  A  total  of  1000  grid  points  in  the  r2 
direction  with  ^r2=0.1  is  chosen.  The  function  and  derivatives  must  remain  zero  at  a 
significant  number  of  nodes  in  order  to  ensure  that  conditions  at  the  computational 
boundary  do  not  constrain  the  solution. 

At  each  value  of  the  “parameter”  t,  the  integration  is  initiated  with  the  initial 
conditions  in  (98),  subject  to  the  boundary  conditions  in  (99)  .  The  spatial  distribu¬ 
tion  of  the  solution  with  respect  to  rj  evolves  as  the  “time-like”  variable  rj  increases. 
Integration  is  carried  out  to  a  sufficiently  large  value  of  rj  to  ensure  that  adequate 

data  fields  Wo{z,t,ri,r2)  are  available  on  the  locus  curve  relating  rj  and  r2  defined 
by  (51).  Then  the  actual  solution  is  found  from  the  intersection  of  the  surface  defined 

by  l^o(2i<)’’ii^j)  and  the  vertical  plane  from  the  locus  curve  relating  ri  and  rj. 

The  first  test  case  studied  is  for  u;  =  1.0,  which  is  a  nonresonant  frequency  smaller 
than  the  first  natural  frequency  Ai  =  |.  This  particular  case  allows  us  to  develop 
a  relatively  simple  solution  with  minimal  computational  time.  The  initial  condition 
was  obtained  from  (98)  by  using  the  first  20  Fourier  modes. 

It  is  important  to  resolve  the  solution  in  the  axial  direction  by  choosing  a  suffi¬ 
ciently  large  number  of  grid  points  in  the  z-direction,  Kmam-  Solution  comparisons  for 
Knuix  =  9|21  and  28  at  t  =  40,  M=0.01  and  Z=0.25,  0.5  and  0.75  show  considerable 
difference  for  the  smallest  number.  Excellent  comparisons  for  Kma*  =  21  and  28 
suggest  that  the  former  is  adequate  for  u;  =  1.0  . 

Figures  6,7  and  8b  give  results  for  the  rotational  axial  velocity  \Vo  ^  zl  function 
of  ri  at  z=0.5  when  t=  20,30  and  40.  At  each  time,  the  solution  consists  of  several 
spatial  oscillations  and  a  thin  region  of  exponential  damping  near  a  specific  radial 
location  ri,ff)  beyond  which  vorticity  is  exponentially  small.  For  a  given  value  of  t, 
one  first  finas  r2,  from  rf  =  0  =  t  —  r2t  and  then  uses  the  inverse  of  (51)  to  determine 
Tig.  It  should  also  be  noted  from  the  definition  of  rj  in  (86)  and  the  second  of  (51)  that 

the  front  moves  at  the  nondimensional  speed  |j|  =  —  MKo*(»’)  which  corresponds 

to  the  steady  radial  speed  in  dimensional  terms.  The  front  speed  vanishes  as  the 
cylinder  axis  is  approached.  At  t  =  40  the  rotational  flow  distribution  in  Figure  4b 
extends  out  about  44%  of  the  cylinder  radius. 

The  corresponding  scaled  vorticity  distribution,  calculated  from  is  given 

in  Figure  9.  The  dimensionless  vorticity  fl*,  defined  below  (84)  is  0(M"^).  One 
should  note  the  significant  magnitude  of  the  spatial  variations  in  vorticity  given  the 
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scaling  factor  The  large  amplitude  is  associated  with  the  0(M)  length  scale  of 
the  axial  speed  gradient  and  the  factor  A/“’ 

The  spatial  distribution  of  the  axial  rotational  velocity  at  each  time  can  be  ex¬ 
plained  in  physical  terms  by  considering  the  interaction  between  the  propagating 
planar  acoustic  wave  solution,  described  by  (69)  and  (70),  and  the  steady  injected 
flow  field  in  (15).  Fluid  particles  injected  normally  from  the  wall  with  =  1  at 
a  specified  z-location  experience  an  approximately  harmonic  variation  in  the  local 
pressure  gradient  given  by  (70).  During  periods  of  negative  (positive)gradients  the 
particles  are  accelerated  downstream  (upstream).  As  a  result,  near  the  w-ali  one  will 
observe  alternating  periods  of  positive  and  negative  axial  velocity.  The  steady  radial 
velocity  field  carries  these  alternating  regions  of  forward  and  reverse  flow  away  from 
the  wall  toward  the  axis.  Part  of  the  fluid  particle  response  is  purely  acoustic.  The 
remainder  is  given  by  Wq.  In  this  sense  the  spatial  pattern  of  the  rotational  axisd  ve¬ 
locity  at  fixed  z  as  shown  in  Figures  6,7  and  8b  reflects  the  histoncal  behaviour  of  the 
local  pressure  gradient  at  the  axial  location.  The  0(M)  length  scale  of  the  transverse 
spatial  oscillations  can  now  be  explained  easily  since  the  approximately  harmonic 
pressure  gradient  variation  occurs  on  the  acoustic  time  scale  (for  u^=0(l))  during 
which  only  limited  radial  motion  is  possible. 

Three  additional  results  should  be  noted.  First,  the  amplitude  of  a  given  oscillation 
in  the  vortical  velocity  component  becomes  smaller  as  is  convected  deeper  into  the 
cylinder.  Local  viscous  and  nonlinear  effects  have  some  influence  on  the  amplitude 
reduction  which  must  occur  because  the  vorticity  vanishes  near  the  front  and  certainly 
at  the  axis  ri  =  1.  Second,  the  wave  length  of  the  spatial  oscillations  decrease  as 
Ti  increases  toward  the  axis.  This  occurs  because  the  vorticity  front  speed  (the 
steady  radial  velocity)  declines  as  ri  increases.  Finally,  a  comparison  of  Figures 
8a,b,c  demonstrates  that  the  rotaion^  axial  velocity  varies  with  the  axial  as  well  as 
the  radial  variables. 

Figures  lOa-c  show  a  “linear”  solution  for  the  rotational  axial  velocity  Wq  as  a 
function  of  ri  at  t=40,  M=0.01  and  z=0.25,  0.5  and  0.75  obtained  by  solving  (97) 

with  nonlinear  term  reduced  by  a  factor  of  10“®.  K^ax  =  21  is  used  in  the 

calculation.  This  calculation  is  done  to  assess  the  impact  of  the  nonlinear  convection 
term  on  the  spatial  structure  of  the  solution  in  Figures  8a-c.  A  comparison  shows 
that  the  nonlinear  term  has  a  quantitative  effect  on  the  spatial  distribution  of  Wq, 
but  does  not  fundamentally  control  the  qualitative  spatial  oscillations.  The  nonlinear 
effect  is  more  important  in  the  fore  end  at  z  =  0.25  than  in  the  rear  end  at  z  =  0.75 
near  the  the  exit  where  nonlinear  effect  disappears  due  to  the  pressure  node. 

Figures  11a,  c  are  the  counterparts  to  Figures  8a,  c  with  a  reduced  viscous  effect. 

In  this  case  (97)  is  solved  with  the  viscous  term,  multiplied  by  a  factor  of 

•*  1  .  ^  ^7 

0.25.  Tne  basic  patterns  of  the  complete  solution  in  Figure  8  persist.  One  notes  that 
oscillation  peak  amplitude  in  Figure  7a,c  are  considerably  larger  than  their  analogues 
in  Figure  8a,  c,  particularly  away  from  the  wall.  Smaller  differences  are  seen  in 
Figure  7b.  In  general,  the  impact  of  viscosity  is  greater  at  an  axied  location  where 
the  maximum  oscillation  amplitudes  are  relatively  large. 

Figure  12  and  Figure  13  are  the  counterparts  of  the  nonlinear  result  in  Figure 
8b  with  larger  axial  Mach  numbers  M  -0.05  and  0.1  respectively,  and  t=40.  One 
can  see  that  the  larger  Mach  numbers  (stronger  wedl  injection)  are  associated  with 
a  thicker  unsteady  vortical  region.  The  vorticity  has  filled  the  entire  cylinder  for 
M=0.05  and  for  M=0.1  at  t=40  while  only  44%  of  the  chamber  is  filled  in  Figure 
8  with  M=0.01.  This  is  expected  because  radial  convection  of  vorticity  occurs  more 
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quickly  for  a  relatively  larger  steady  radial  velocity.  Although  vorticity  is  present 
everywhere  in  the  cylinder,  the  value  at  the  axis  is  zero,  as  explained  below  (82).  It 
should  be  noted  that  the  local  velocity  gradients  are  smaller  in  Figures  12  and  13 
so  that  the  magnitude  of  the  unsteady  vorticity  is  similarly  reduced  for  higher  Mach 
number  systems.  Larger  velocity  gradients  can  be  obtained  at  a  given  Mach  number 
for  larger  forcing  frequency  w. 

Figure  14  provides  results  at  ri  =  0.2  for  the  time-variation  of  the  planar  acoustic 
audeJ  speed  VFopi  the  rotational  axial  speed  Wq  and  their  sum  V^o,  defined  in  (59),  at 
z=0.5  with  Kmax  =  21  and  M=0.01.  It  can  be  seen  from  Figure  14b  that  the  curve 
for  the  rotational  solution  resembles  the  planar  acoustic  response  in  Figure  14a  but 
differs  in  amplitude  and  phase.  At  ti  =  0.2  the  rotationzd  response  appears  after  a 
delay  of  almost  18  axial  acoustic  time  units,  the  time  needed  for  the  vorticity  wave 
front  initiated  at  the  wall  to  travel  out  to  the  specified  radial  location.  At  the  location 
ri  =0.2,  phase  differences  are  relatively  small  and  the  sum  in  Figure  14c  shows  a  total 
response  of  significant  amplitude.  This  amplitude  actually  increases  as  ri  decreases 
until  locations  very  close  to  the  wall  are  reached,  where  the  impact  of  the  no-slip 
condition  at  ri  =  rj  =  0  forces  V^o  — *  0. 

The  second  case  studied  is  for  a  driving  frequency  u  =  2.5  with  the  initial  condition 
for  Wo  in  (98)  half  the  magnitude  of  the  previous  calculation.  The  solution  for  Wo 
at  t=40,  M=0.01  is  given  in  Figures  15a-c  for  ifmoa  =  21  .  The  results  show  that  the 
spatial  distribution  curves  for  a;  =  1.0  and  2.5  have  similar  cheiracteristics.  However, 
the  latter  have  more  oscillation  cycles  because  the  driving  frequency  is  higher.  One 

should  note  that  the  axial  variation  of  Wq  is  quite  different  from  the  lower  frequency 
analogue  in  Figure  8a-c. 

The  relative  solution  complexity  for  u;  =  2.5  suggests  that  a  less  restricted  acoustic 
field,  arising  from  multiple  driving  frequencies  or  perhaps  sidewall  injection  oscilla¬ 
tions,  may  initiate  a  relatively  irregular  vortical  time-response.  In  this  sense,  one 
could  aisk  whether  “turbulent”  responses  observed  in  similar  situations  such  as  solid 
rocket  chamber  models  aire  caused  in  part  by  wall  generated  vorticity  that  is  convected 
into  the  chamber  by  the  injected  field. 


7  Modal  Solutions  to  the  Nonlinear  System 

Solutions  found  from  (101)-(107)  enable  one  to  demonstrate  how  the  energy  is  dis¬ 
tributed  among  the  spectrum  of  Fourier  modes.  The  coefficients  An{ri,r2),  n=l-N 
for  specified  N,  have  been  found  by  using  a  finite  difference  method  based  on  the 
Adam-Bashforth/Crank-Nicolson  scheme  described  in  Section  6.  The  procedure  used 

previously  to  find  Wq  oil  the  real  locus  is  applied  to  each  A„(ri, r2).  Once  the  A'^a  are 

known,  the  solution  for  Wo  is  found  by  summing  up  N  modal  contributions  according 
to  (101)';  Careful  attention  must  be  given  to  the  value  of  N  in  order  to  assure  that 

the  solutions  are  sufficiently  accurate.  In  particular,  solutions  for  Wo  based  on  N  and 
N  +  L  modes,  L  >  0,  are  compared  until  the  results  exhibit  little  appreciable  change 
for  an  incremental  L  value.  V2duable  insights  into  the  solution  development  for  this 
problem  have  been  found  from  the  transient  solution  to  a  Fisher  equation  (Fisher, 
1936)  in  terms  of  a  Galerkin  expansion.  Details  are  given  in  Appendix  A. 

Figures  16-18  show  Wq  results  at  t=40  for  N  =  6,8,10  respectively  for  the  case 
M=0.01,z=0.5  and  u;  =  1.0.  Each  figure  includes  partial  sums  of  even  numbers  of 
modes  up  to  the  complete  finite  summation.  A  comparison  of  partial  sum  results 
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implies  that  the  energy  is  concentrated  primarily  in  the  first  two  modes  for  w  =  1.  A 
comparison  of  results  in  Figures  16-18  with  those  in  Figure  8b  reveals  that  there  is 
little  difference  between  the  modal  solutions  and  that  from  a  direct  finite  difference 
solution  with  K^ax  =  21.  It  is  noted  that  the  grid  size  in  the  axial  direction  for 
l^mo*=21  is  0.05,  just  small  enough  to  resolve  the  first  six  axiad  Fourier  modes  ac¬ 
cording  to  {Az)max  =  where  Aa  =  ~ir.  Hence  the  most  meaningful  comparisons 
should  be  carried  out  with  iV  =  6. 

It  has  been  found  that  the  6  mode  partiad  sum  from  the  10  mode  computation 
is  closer  to  the  finite  difference  result  thaji  the  6  mode  summation  from  the  6  mode 
computation.  A  rea.sonable  explanation  of  this  observation  may  be  baised  on  the 
restricted  energy  transfer  for  modes  neaj  the  truncation  limit.  Adding  a  few  more 
modes  permits  readistic  exchange  between  the  lower  modes.  Hence  a  6  mode  partial 
sum  from  an  Af  =  10  calculation  provides  better  results  thain  the  N  =  6  calculation 
alone.  This  effect  occurs  in  a  related  work  by  Wang  and  Kaissoy  (1993). 

Similar  modal  computations  done  for  u  =  2.5  with  the  initial  mode  amplitude 
reduced  by  50%  are  shown  in  Figure  18  for  N  =  6.  Rewonably  good  comparisons  are 
found  from  a  six  mode  summation.  Partial  sum  comparisons  suggest  that  the  energy 
is  again  concentrated  in  the  first  few  modes. 

The  smaller  amplitude  vortical  radial  velocity  can  be  calculated  by  integrating 
equation  (77).  Figure  19  describes  the  rotationid  radial  velocity  as  a  function  of 
n  at  t=40,  M=0.01  and  z—  0.5  with  a  6  mode  partial  sum  from  a  8  mode  computa¬ 
tion.  Beyond  the  vorticity  front  the  transient  radial  velocity  distribution  is  a  known 
constant  multipl3ring  the  steady  radiad  velocity. 

The  satisfying  comparison  between  finite  difference  and  modal  solutions  provides 
strong  confidence  in  the  characteristic  solution  properties.  It  is  noted  that  CPU 
time  requirement  for  direct  finite  diffence  computations  is  considerably  less  than  that 
for  modal  solutions,  given  equivalent  resolution  requirements.  However,  the  modal 
solutions  can  be  used  to  obtain  insights  to  the  energy  distribution  at  various  length 
scales,  not  easy  to  find  from  finite  difference  calculations. 


8  Conclusions 

Systematic  methods  have  been  employed  to  formulate  a  mathematical  model  to  de¬ 
scribe  the  creation  and  evolution  of  unsteady  rotational  flow  in  a  long,  narrow  cylin¬ 
der.  Boundary  driven  axial,  planar  acoustic  waves  interact  with  an  invisdd,  weakly 
rotational,  injection  induced  steady  flow  to  produce  time  dependent  vorticity  at  the 
sidewall  of  the  cylinder.  The  relatively  intense  vorticity  is  convected  into  the  en¬ 
tire  chamber  by  the  steady  radial  velocity  field  for  appropriate  ranges  of  Reynolds 
and  Mach  number  and  frequency.  The  amplitude  and  distribution  of  the  vorticity  is 
impacted  by  viscous  and  nonlinear  effects. 

It  is  'alio  demonstrated  that  there  are  parameter  ranges  of  Mach  number  (injection 
rate),  driving  frequency  and  Reynolds  number  for  which  vorticity  is  really  confined  to 
thin,  classics,  but  injected  acoustic  boundary  layers  like  those  discussed  by  Flandro 
(1974),  Baum  and  Levine  (1987)  .  These  structures  can  appear  for  relatively  small 
injection  rates,  relatively  high  driving  frequency  and  low  Reynolds  numbers,  so  that 
viscous  damping  of  the  vorticity  amplitude  is  profound.  Then,  the  cylinder  core 
will  contain  the  relatively  weak  vorticity  of  the  steady  Culick  (1966a)  solution  and 
irrotationai  acoustic  waves  driven  by  the  boundary  forcing. 

There  is  now  a  considerable  body  of  evidence  in  support  of  the  presence  of  an 
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unsteady  vorticity  distribution  within  an  appropriately  high  Reynolds  number  wall 
injected  flow  in  a  cylinder.  The  experiments  of  Brown  et  al.(1986a, 1986b),  the  quasi- 
analytical  modeling  of  Flandro  and  Roach  (1993)  and  the  computational  solutions 
of  Vuiilot  and  AvaJon  (1991), Flandro  and  Roach  (1993),  and  that  of  Kirkkopru  et 
al.  (1994)  as  well  as  the  current  work  show  unequivocally  that  unsteady  vorticity 
is  generated  at  the  cylindrical  surface  and  is  convected  away  by  the  injected  fluid. 
The  core  of  the  cylinder  is  free  of  vorticity  only  during  the  very  early  phases  of  the 
transient  process,  prior  to  the  arrival  of  a  well  deflned  unsteady  vorticity  front. 

The  amplitude  of  the  transient  vorticity  distributions  described  by  Kirkkopru  et 
al.(1994),  and  in  the  present  work  are  larger  than  that  of  the  Culick(  1966a) 

steady  solution.  It  follows  that  there  will  be  a  relatively  large  axial  shear  stress  on  the 
cylinder  surface, which  can  be  calculated  from  equation  (84),  particularly  for  smaller 
M  values.  This  result  is  important  for  applications  of  the  theory  to  solid  rocket 
engines. 

One  can  speculate  that  the  large  transient  shear  stresses  will  impact  the  burning 
rate  of  a  propellant  which  is  the  source  of  the  “injected"  fluid  used  in  the  present 
model.  Perhaps  there  is  a  direct  relationship  between  the  effect  of  surface  shear  stress 
transients,  predicted  in  the  present  work,  and  erosive  burning  concepts  used  in  the 
solid  rocket  engineering  literature  (Williams, 1985). 

One  may  reasonably  conclude  from  the  results  of  the  aforementioned  studies  that 
the  classicsd  acoustic  stability  predictions  used  for  scientiflc  and  engineering  pur- 
poses(e.g.  Culick,1990)  should  be  reassessed  to  determine  the  prediction  reliability. 
It  is  conceivable  that  the  acoustic  pressure  predictions  are  reasonable  because  the 
presence  of  vorticity  does  not  impact  the  primary  acoustic  pressure  transient  found 
in  the  current  theory.  However,  it  is  unlikely  that  velocity  and  shear  stress  predictions 
are  accurate,  given  the  irrotational  basis  of  the  acoustic  stability  analysis. 

The  conceptual  approach  used  here  can  be  extended  to  disturbances  driven  by 
sidewall  injection  transients,  rather  than  those  applied  at  the  closed  endwall.  The 
former  type  of  disturbance  emulates  the  effects  of  propellant  burning  rate  variations 
in  solid  rocket  engines.  Additionally  it  is  noted  that  the  variable  scaling  reported  here 
will  enable  the  development  of  more  accurate  numerical  methods  like  those  described 
in  Kirkkopru  et  al.(1994). 
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Appendix  A:  Solution  to  a  Related  Model  Problem 


The  nonlinear  coupled  system  in  (102)  is  sufficiently  complex  to  require  a  computa¬ 
tional  solution.  In  order  to  develop  an  effective  numerical  approach,  it  is  desirable 
to  consider  the  solution  to  an  elementary  model  problem  with  related  properties.  A 
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simple  Fisher  equation  (Fisher, 1936)  with  appropriate  periodic  initial  and  boundary 
conditions  can  be  used: 


ds 


dy^ 


+  uU\ 


Initial  Condition; 


y  >  0,t  >  0 


y)  =  -  sin(t  -  y)  for  0  <  y  <  t; 
^(0,y)  =  0  for  y  >  t 

Boundary  Condition: 

U{y  =  0)  =  -  sin(t) 


{A-l) 

{A -2) 
{A -3) 

(A -4) 


When  the  parameter  t  is  increased,  the  nonzero  portion  of  the  initial  condition  is 
spread  farther  into  the  y-domain.  In  the  spirit  of  (102),  the  multiple  scale  independent 
variables  are  related  by  y  =  and  fl  1. 

An  analytical  solution  for  linear  diffusion  (t/  =  0)  is  constructed  for  the  odd  ex¬ 
tension  of  (A-2)  and  (A-3)  for  the  domain  0  <  y  <  oo: 


Uis,y) 


f  sin(t  -  y') 
2-)/irs  Jo 

ysin(f)  r* 

2y/v  Jo  (j  -  a»)l 


e  —  e 

ds',  y  >  0 


dy' 


(A -5) 


A  quasi-steady  solution  form  U{s,y)  —  —  sin(t  —  y)c"*,  y  >  0  can  be  recovered  by 
taking  the  limits  (jf^)  — ►  ±oo  and  s  — »  O'*'  simultanously.  Physically,  this  means 

that  -the  solution  has  a  quasi-steady  form  at  a  specific  value  of  t  if  y  lies  between  y=0 
and  the  inner  edge  of  a  diffusive  boundary  layer  centered  at  y  =  t  which  is  ne^ed 
to  smooth  the  discontinuous  slope  of  the  initial  condition  (A-2)  and  (A-3)  at  that 
location.  Inside  the  diffusive  layer,  the  solution  is  given  by  the  full  form  of  (A-5). 

The  diffusive  layer  thickness  is  ^  ^(^)- 

The  solution  to  (A-l)-(A-4)  along  the  locus  y  =  Q^s,  for  */  =  1  and  =s  66, has 
been  found  from  a  computational  analysis  based  on  an  elementary  explicit  finite 
different  method.  The  boundary  condition  (A-4)  is  enforced  at  y  =  0  for  each  in¬ 
tegration  step  in  s  direction.  The  integration  of  [/  with  respect  to  a  is  carried  out 
for  0  <  a  <  2.25  with  step  size  Sa  =  0.0015.  The  far  end  boundary  condition  is 
implemented  in  such  a  way  that  there  are  a  sufficient  number  of  grid  points  with  zero 
value  lying  between  the  furthest  grid  point  with  nonzero  value  and  the  finite  location 
of  the  computational  boundary  with  respect  to  y  after  each  integration  in  a  direction. 
The  dashed  line  in  Figure  20  describes  the  linear  solution  when  t  =  100  obtained  from 
(A-5).  An  analogous  numerical  result  (i/  =  0)  is  indistinguishable  from  the  analytical 
solution  on  the  scale  of  the  graph,  thus  verifying  the  numerical  code.  The  linear  solu¬ 
tion  shows  regular,  nearly  harmonic  spatial  oscillations  that  decay  until  the  diffusive 
layer  is  reached  near  a  %  1.5.  There  the  solution  makes  a  rapid  transition  to  a  van- 
isUngly  small  value  for  s  >  1.5.  In  comparison  the  solid  line  represents  the  nonlinear 
numerical  solution  for  t/  =  1.  The  frequency  is  nearly  identical  to  the  linear  solution. 
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However,  the  drift  of  the  solution  toward  positive  values  of  U  is  due  to  the  positive 
definite  source  effect,  i/U^.  Again  the  deviation  from  the  pattern  of  oscillations  near 
3  =5:  1.5  is  associated  with  the  diffusive  layer  behavior.  Given  the  parameters  used  in 
the  calculation,  the  diffusive  layer  thickness  with  respect  to  the  s  coordinate  is  about 
0.1.  The  analogous  results  for  t/  =  —  1,  corresponding  to  a  nonlinear  sink,  are  given 
in  Figure  21.  There  is  no  expectation  of  symmetry. 

The  basic  properties  of  the  model  problem  solution  can  be  used  to  develop  an 
effective  numerical  method  for  (102). 
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Table  1:  Acoustic  Response  Properties  for  Several 
Driving  Frequencies 


u; 

u;'(Hz) 

Properties 

Primary  Response 

1 

159 

stable 

axial  -1-  quasi-steady  modes 

238 

beats 

quasi-steady  modes 
with  axial  wave  modulation 

ir/2 

250 

axial  amplification 

linear  growth 
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9  CAPTIONS 


Figure  1:  The  cylindrical  rocket  engime  chamber  model  of  length  L' ,  diameter  D' , 
with  endwall  oscillations  of  frequency  u>. 


Figure  2:  The  time  response  of  the  ajdal  acoustic  velocity 


Figure  3:  The  time  response  (beats)  of  the  axial  acoustic  velocity  K,o/30  at  z=0.5 
and  a;  =  1.5. 


Figure  4:  The  time  response  of  the  axial  acoustic  velocity  at  z=0.5  and  =■ 
ir/2, linear  growth. 


Figure  5:  The  axial  velocity  in  the  boundary  layer  as  a  function  of  the  variable  of 
(  fo  M=0.01,/3  =  0.1  when  fi  =  2.5  and  3.0  .  The  boundary  layer  is  thicker  for  the 
smaller  frequency  value. 


Figure  6:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wo  with 
the  radial  variable  vi  at  t=20.  It  can  be  seen  that  the  front  is  about  23%  of  the  way 
to  the  centerline  and  nearly  4  spatial  oscillations  have  entered  the  cylinder. 


Figure  7a,b,c:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wo  at 
t=30  for  (a)z=0.25,(b)z=0.5  fc)z=0.75  with  the  radial  variableri.  The  front  is  about 
35%  of  the  way  to  the  centerline  and  nearly  5.5  spatial  oscillations  have  entered  the 
cylinder.  The  amplitude  of  the  oscillations  near  the  front  at  this  moment  is  noticably 
smaller  those  at  t=20,  an  accumulative  effect  of  viscosity  on  M  scale. 


Figure  8a:  The  spatial  variation  of  the  rotational  axial  velocity  component  Vt’o  "dth 
the  radial  vsnriable  rj  at  t=40,  z=0.25.  The  front  is  about  45%  of  the  way  to  the  cen¬ 
terline  and  nearly  7.5  spatial  oscillations  have  entered  the  cylinder.  Strong  nonlinear 
effects  alter  the  overall  pattern. 
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Figure  8b:  The  spatial  variation  of  the  rotational  axial  velocity  component  Hg  with 
the  radial  variable  rj  at  t=40,  z=0.5.  The  front  is  about  45%  of  the  way  to  the 
centerline  and  nearly  7.5  spatial  oscillations  have  entered  the  cylinder.  The  amplitude 
of  the  oscillations  near  the  front  at  this  moment  is  substantially  smaller  those  at  t=20, 
an  accumulative  effect  of  viscosity  on  M  scale. 


Figure  8c:  The  spatial  variation  of  the  rotational  axial  velocity  component  VVo  with 
the  radial  variable  ri  at  t=40,  z=0.75  .  It  has  the  generic  features  as  in  Figure  8b. 


Figure  9:  The  spatial  variation  of  unsteady  vorticity  fij/lOO,  with  the  radial  variable 
rj  .  The  magnitude  of  Qs  is  0(1/M)  bigger  than  the  steady  vorticity  and  therefore 
dominant  in  the  cylinder.  The  shear  stress  on  the  wall  imposed  by  the  unsteady 
vorticity  is  significant. 


Figure  10a, b,c:  The  spatial  variation  of  the  rotational  axiad  velocity  component  Wq  at 
t=40  for  (a)z=0.25,(b)z=0.5,(c)z=0.75,  with  the  radial  variable  rj  when  the  nonlinear 
term  is  suppressed.  A  comparison  of  these  with  Figures  8a,b,c  demonstrates  nonlinear 
effects  are  strong  near  the  fore  end,  moderate  at  the  mid^e,  small  at  the  rear  end. 


Figure  lla,b,c:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wo  at 
t=40  for  (a)z=0.25,^)2=0.5,(c)z=0.75,  with  the  radial  variable  rj  with  the  viscous 
term  reducM  by  50%.  A  comparison  with  Figures  8a,b,c  demonstrates  that  local 
structure  is  altered  by  the  reduction  of  viscosity. 


Figure  12:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wo  with 
the  radial  variable  tx  at  M=0.05,  z=0.5,  u;=1.0  and  t=40.  The  larger  M  corresponds 
to  stronger  injection  rate  on  the  cylinder  wall.  As  a  result,  almost  all  the  cylinder 
has  been  filled  with  the  rotational  flow  compared  to  only  45%  in  the  case  of  M=0.01 
in  Figure  8c. 


Figure  13:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wo  with 
the  radial  variable  rx  at  M=0.1,  z=0.5,  u;=1.0  and  t=40.  The  larger  M  corresponds 
to  a  stronger  injection  rate  on  the  cylinder  wall. 


Figure  14  a,b,c  (&om  top  to  buttom):  The  time  response  of;  (a)  the  axial  acoustic 
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speed  VVop,(b)  the  rotational  ajcial  speed  U'o  and  (c)  the  complete  axial  speed  V^o  = 
(IVop  +  Wo)  for  ri=0.2, 2=0.5,  M=0.01  and  u;  =  1.0. 


Figure  15a, b,c;  The  spatial  variation  of  the  rotational  axiad  velocity  component  W'o 
with  the  radial  variable  rj  at  u;  =  2.5,t=40  for  {a)z=0.25,(b)z=0.5,  (c)z=0.75  with 
the  amplitude  of  the  disturbance  reduced  by  a  half.  More  spatial  oscillations  are 
present  in  the  structure  due  to  the  higher  driving  frequency,  in  comparison  Figures 
8a,b,c. 


Figure  16:  The  spatial  variation  of  the  rotational  aixial  velocity  component  Wq  with 
the  radial  variable  ri  based  on  a  6  mode  summation  from  a  6  mode,  nonhnear  com¬ 
putation  (N=6). 


Figure  17:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wq  with 
the  radial  variable  ri  based  on  a  6  mode  partial  summation  from  an  8  mode,  nonlinear 
computation  (N=8). 


Figure  18:  The  spatial  variation  of  the  rotational  axial  velocity  component  Wo  with 
the  radial  variable  ri  based  on  a  6  mode  partial  summation  from  a  10  mode,  nonlinear 
computation  (N=10). 


Figure  19:  The  spatial  variation  of  the  rotational  radial  velocity  component  with 
the  radial  variable  ri  based  on  the  6  mode  partial  summation  from  am  8  mode  com¬ 
putation. 


Figure  20:  Solution  U  vs  z  for  the  nonlinear  model  problem  on  the  locus  curve  y  = 
with  n*  J5S  66  and  v  =  1  The  dashed  line  is  the  plot  of  U  vs  z  from  (A-5)  on  the  same 
locus  curve. 


Figure  21:  Solution  U  vs  z  for  the  nonlinear  model  problem  on  the  locus  curve  y  =  Cl^z 
with  fl*  %  66  and  v  =  —1  The  dawhed  line  is  the  plot  of  U  vs  z  from  (A-5)  on  the 
same  locus  curve. 
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Figure  1 :  The  cylindrical  rocket  engine  chamber 
model  of  length  L' ,  diameter  D' ,  with  endwall 
oscillations  offt'equency  o) 


Appendix  c 


Unsteady  Vorticity  Generation  and  E\olution  in  a 
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ABSTRACT 


Two  dimensional,  ajcisymmetric  Navier-Stokes  equations  are  solved  numerically  to 
study  the  unsteady,  nonlinear  vorticaJ  field  generation  and  evolution  by  an  acoustic 
wave  interaction  with  a  sidewall  injected  flow  in  a  cylindricaJ  tube,  that  mimics  the 
solid  propellaint  surface  burning  in  a  rocket  engine.  The  steady  flow  field,  sustained 
by  the  sidewall  injection,  is  perturbed  by  imposing  a  sinusoidal  component  on  the  exit 
plane  static  pressure.  Amplitudes  of  the  oscillatory  pressure  disturbances  are  chosen 
accordingly  for  different  injection/  mean  axial  flow  Mach  numbers  so  that  nonlinear 
processes  affect  the  evolution  of  the  vortical  field.  In  the  spirit  of  the  study  by  Zhao 
et  al.*®,  the  unsteady  vortical  part  of  the  total  velocity  field  is  extracted  from  the 
numerical  solutions  to  show  explicitly  the  vorticity  generation  at  the  injecting  side 
wall  and  evolution  in  time.  Present  results  show  the  unsteady  vorticity  eventuaDy 
fills  the  entire  chamber  unlike  that  confined  into  traditional  acoustic  viscous  layers 
adjacent  to  the  sidewall. 

I.  INTRODUCTION 

Transient  flow  dynamics  in  a  solid  rocket  engine  chamber  are  strongly  coupled  to 
propellant  combustion  processes  and  related  directly  to  overall  rocket  motor  stability. 
Traditionally,  prediction  models  have  been  formulated  in  terms  of  acoustic  stability 
theory.  The  earliest  study  was  devised  by  Grad^.  Subsequently,  Culick  and  co- 
workers^  developed  linear  and  weakly  nonlinear  stability  theories  to  describe  the 
response  of  model  systems  to  assumed  acoustic  disturbances.  At  about  the  same  time. 
Hart  and  McClure^^  proposed  a  theory  based  on  the  acoustic  balance  approach  for 
prediction  of  the  flow  stability  in  the  rocket  engine  chamber.  This  theory,  based  on 
linearized  invisdd  equations,  has  been  employed  to  evaluate  laboratory  and  full-size 
rocket  engine  performance.  Prediction  of  pressure  responses  compares  well  with  those 
measured  in  experiments. 

Recently,  Brown  et  conducted  laboratory  experiments  in  a  cold  flow  rocket 
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engine  chamber  analogue.  Velocity  measurements  taken  along  and  across  the  cylin¬ 
drical  chamber  show  that  there  is  a  significant  unsteady  rotational  flow  component 
present  everywhere.  This  type  of  vorticity  is  seen  adso  in  the  study  by  Vuillot  and 
Avalon*®  who  used  computational  methods  to  solve  the  compressible  Navier-Stokes 
equations  in  a  channel  with  side  wall  mass  injection.  These  experimental  and  numeri¬ 
cal  results  suggest  that  stability  predictions  based  on  the  traditional  acoustic  analysis 
should  be  reexamined. 

Traditionally,  stability  predictions  are  based  on  mathematical  models  formulated 
in  terms  of  linear  or  weakly  nonlinear,  invisdd,  irrotational  acoustic  concepts. 
(  Culick*"**, Williams*®).  In  these  models,  acoustic  waves  propagate  through  a  quies¬ 
cent  chamber,  and  do  not  interact  with  side  wall  gas  injection  that  mimics  maas  input 
to  the  chamber  by  solid  propellant  burning.  These  mathematical  models  cannot  de¬ 
scribe  the  vorticity  distributions  observed  in  laboratory  and  numerical  experiments, 
which  are  discussed  in  a  Umited  way  by  investigators  concerned  with  classical  acoustic 
boundary  layers.  (  for  example,  Flandro*^,  Wang  and  Kassoy*^). 

Flandro*^  demonstrated  that  rotational  flow  effects  are  of  considerable  importance 
in  acoustic  boundary  layers  where  viscosity  is  significzmt.  Vmllot  zmd  Kuentzman*® 
employed  Flandro’s  analysis  to  evaluate  a  laboratory  scale  rocket  engine  ,  and  ob¬ 
tained  good  comparisons  between  theory  and  experiment.  The  main  objective  of  these 
studies  was  to  understand  energy  exchange  between  acoustic  disturbances  and  the 
mean  flow  as  the  normally  injected  fluid  is  turned  toward  the  axial  direction.  Baum 
and  Levine**  and  Baum*®  numerically  solved  the  complete  Navier-Stokes  equations  in 
order  to  study  energy  exchange  mechanisms  in  thin,  viscous  acoustic  boundary  layers. 
The  numerical  results  of  Vuillot  and  Avalon*®  demonstrate  that  unsteady  vorticity 
is  not  eilways  confined  to  thin  viscous  acoustic  boundary  layers  adjacent  to  the  in¬ 
jecting  wall.  Instead,  the  unsteady  vorticity  generated  at  the  wall  is  convected  away 
from  the  wall  into  the  chamber,  and  persists  over  a  significant  part  of  the  chamber 
for  appropriate  parameter  values.  The  presence  of  vorticity  in  the  chamber  must  be 
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accounted  for  in  an  accurate  and  reliable  theory  of  engine  stability. 

Flandro  and  Roach^^  developed  an  approximate  model  and  an  analytical  solution 
to  describe  vorticity  generation  at  the  injecting  wall  due  to  an  interaction  of  axial 
acoustic  waves  and  the  injected  flow.  The  model  is  baised  on  purely  invisdd  equations. 
An  effort  was  also  made  to  simulate  numerically  the  Brown  et  al.^^'^^  laboratory 
experiments  by  solving  the  Navier-Stokes  equations.  This  work  implies  that  there 
are  two  important  length  scales  for  rotational  effects;  the  tube  radius  and  an  0{M) 
smaller  length  where  significant  local  velocity  variations  are  present. 

Recently,  Zhao  and  Kassoy**  as  well  as  Zhao  et  al.**  used  systematic  asymptotic 
methods  to  formulate  an  initial  boundary  value  problem  that  includes  rotational  flow 
effects.  In  particular,  the  theory  describes  the  generation  amd  evolution  of  unsteady 
vorticity  produced  at  the  cylinder  wall  by  an  interaction  between  the  injected  fluid 
and  the  propagating  planar  acoustic  disturbances.  It  is  demonstrated  analytically, 
to  a  first  order  approximation  for  a  large  Re  and  low  M  flow,  that  the  interaction 
between  the  acoustic  disturbance  and  the  injecting  fluid  is  purely  in\iscid.  Analysis 
is  also  used  to  show  that  the  vorticity  generated  at  the  wall  is  convected  out  into 
the  primarily  invisdd  core  flow  by  the  radial  component  of  the  injection-induced  flow 
field.  This  asymptotic  model  formally  incorporates  two  distinct  length  scales,  similar 
to  those  described  in  the  study  by  Flandro  and  Roach^^,  in  a  multiple  scale  analysis 
which  indudes  weak  viscous  and  nonlinear  effects. 

The  objective  of  the  present  study  is  to  compute  unsteady  vortidty  production  and 
evolutionin  a  finite  cylinder  with  steady  wall  injection  and  an  imposed  disturbance  on 
a  downstream  exit  plane.  Having  invaluable  information  about  the  existence  of  two 
disparate  length  scales  and  their  relationship  to  Mach  number,  enables  grid  size  and 
spatial  distribution  to  be  chosen  carefully  in  order  to  describe  the  flow  dynaunics  accu¬ 
rately.  Compressible  Navier-Stokes  equations  are  solved  by  utilizing  the  MacCormack 
explidt,  predictor-corrector  method”.  A  low  Mach  number,  M  =  0(10“*  —  10“^), 
high  Reynolds  number  Re  =  0(10^  —  10®),  steady  gas  flow  arising  from  constant  side- 
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wall  injection  in  a  long,  narrow  cylindrical  chamber  is  disturbed  at  the  exit  plane  by  a 
sinusoidally  fluctuating  planar  pressure  variation.  Following  the  approach  of  Flandro 
and  Roach*^  Zhao  and  Kassoy*^  and  Zhao  et  the  rotational  component  of  the 
velocity  field  is  extracted  from  the  numerical  solutions  for  the  total  velocity.  The 
computational  results  show  that  unsteady  vorticity  generated  at  the  wall  by  an  in¬ 
teraction  between  -  he  injected  fluid  and  planar  propagating  acoustic  pressure  waves 
is  convected  away  from  the  wall  towards  the  mainly  invisdd  core  flow’  by  the  mean 
radird  flow  velocity  field.  Eventually,  unsteady  vorticity  fills  the  entire  chamber.  The 
presence  of  rotational  flow  features  imply  that  traditional  acoustic  balance  theories, 
used  wridely  to  predict  solid  rocket  engine  chamber  stability,  must  be  reevaluated. 


II.  COMPUTATIONAL  MODEL 


A  primary  goal  of  the  present  numerical  effort  is  to  complement  the  analytical 
study  by  Zhao  et  al.^*  describing  the  unsteady  nonlinear  vorticity  generation  and 
evolution  by  acoustic  wave  interaction  in  a  model  of  a  solid  propellant  rocket  en¬ 
gine  chamber.  For  computational  convenience  a  cylindrical  geometry  is  chosen.  The 
flow  inside  the  chamber  is  sustained  by  constant  velocity  fluid  injection  through  the 
sidewall  of  long,  narrow  cylindrical  tube  with  one  end  closed. 

The  flow  field  is  described  by  axisymmetric,  two-dimensional,  laminar,  compress¬ 
ible  Navier-Stokes  equations  that  are  written  in  nondimeasionail  conservative  form  as 
follows; 
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The  equation  of  state  for  a  perfect  gas  is 

V~pT  (3) 

Nondimensional  variables,  defined  in  terms  of  dimensional  quantities  denoted  by  a 


prime,  are  given  by 
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T  =  nu 

t  =  t'/t'. 

c.  = 

C'JC'^ 

(4) 

Characteristic  length  scales 

for  the  axial  and 

radial  directions  are  chosen  to  be  the 

length  of  the  tube  V  and  the  radius  of  the  tube  i?',  respectively.  The  characteristic 
sidewzdl  injection  speed  of  the  fluid  V|^  is  related  to  the  characteristic  mean  axial  speed 
U'n  through  the  global  mass  conservation  relationship  U'n  =  SVJi  where  6  —  L'/R'  is 
the  aispect  ratio  of  the  tube.  Pressure  is  nondimensionalized  with  respect  to  the 
static  pressure  at  the  open  end  (outlet)  of  the  cylindrical  chamber.  Characteristic 
values  for  temperature  and  density,  and  p'^  respectively,  represent  the  injected 
fluid  properties.  Time  is  nondimensionalized  with  respect  to  the  tube  axial  acoustic 
time  =  L'/uq  where  a'^  =  (7Po/Po)'^*  characteristic  speed  of  sound.  Here, 
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the  ratio  of  specific  heats  7  =  1.4  is  used  in  the  present  computations.  The  viscosity 
specific  heats  and  conductivity  are  treated  as  constants  in  this  calculations  because 
temperature  variations  are  very  small. 

The  following  expressions 


Re 


K  °-'o 


(5) 


define  the  Reynolds  number,  the  Prandtl  number  and  the  mean  axial  flow  Mach 
number,  respectively.  In  a  typical  solid  rocket  engine  chamber  Re  >>  1,  Pr  =  0(1) 
and  M  =  0(10-2  -  IQ-^). 

Finally, 

E,  =  pCJ  +  7(7  -  l)M2l - (6) 

represents  the  nondimensional  form  of  the  total  energy  of  the  fluid. 

The  Navier-Stokes  equations  are  simplified  by  ignoring  the  axial  transport  terms. 
Justification  for  the  reduction  is  based  on  the  asymptotic  anzdysis  in  Zhao  and 
Kassoy22  and  Zhao  et  2l1.2®  valid  for  large  aspect  ratio  5  >>  1  and  the  large  Reynolds 
number  iZe  >>  1  provided  that  P IRe  «  1.  The  computation  time  for  the  simphfied 
equations  is  reduced  drastically  without  sacrificing  the  flow  physics.  Furthermore,  the 
dissipative  effects  of  the  remaining  transport  terms  are  sufficient  to  avoid  artificial 
damping  terms  needed  in  other  similar  computations2^'*^. 

The  Navier-Stokes  equations  are  solved  by  using  the  MacCormack  explicit,  predictor- 
corrector  scheme  : 


At ,  „  „  ,  At ,  „  _  ,  At ,  _ 

9t,j  ~  9i,i  ~  ~  ~  ~ 


Ar' 


1  r 

2 


lij  ~  (^’j  ~  [Jij  ~  f ij-lj  ~  j. 


(7) 


Here,  the  overbar  denotes  the  predictor  stage,  while  the  superscripts  n  and  n  +  1 
represent  the  known  and  unknown  time  levels,  respectively,  separated  by  At.  The 
subscripts  i  and  j  refer  to  axial  and  radial  directions,  respectively. 
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The  upper  half  of  the  cylindrical  chamber  is  taken  as  the  computational  domain 
because  the  flow  is  assumed  axisymmetric.(see  Fig.  1)  A  constant  speed  fluid  is 
injected  into  the  cylinder  through  the  upper  wall.  The  left  side  of  the  chamber  is  a 
closed  rigid  wall  (head  end)  and  the  right  side  of  the  chamber  is  a  flow  exit  where  a 
speciflc  transient  pressure  variation  is  assumed  to  exist. 

The  aspect  ratio  =  20  for  the  present  computations.  Grid  points  are  equally 
spaced  in  each  direction.  Radial  grid  size  is  dependent  on  the  value  of  M  for  reasons 
that  win  be  discussed  in  Section  III. 

Il.a.  Steady  State  Computations 

A  steady  state  flow  solution  is  required  as  an  initial  condition  for  the  transient  flow 
computation.  Boundary  conditions  include  an  impermeable  wall  at  i  =  0  (u  =  0), 
a  static  pressure  condition  at  the  open  end  s  =  1  (p  =  1),  a  specified  injection 
velocity  (u  =  —  1),  temperature  (T  =  1)  and  no  slip  condition  for  the  axial  flow  speed 
(u  =  0)  on  the  injecting  upper  sidewall  at  r  =  1  and  symmetry  conditions  on  the 
lower  (centerline)  boundary,  r  =  0. 

Two  approaches  have  been  used  to  compute  the  steady  state  flow  configuration. 
In  the  first,  the  sidewall  injection  is  initiated  at  t  =  0,  and  the  internal  flow  is  started. 
Numerical  integration  is  carried  out  until  a  specified  convergence  criterium  aissures 
a  steady  state.  For  example,  the  injected  total  mrss  flow  rate  must  equal  the  total 
exiting  mass  flow  rate  and  the  sum  of  changes  in  axial  and  radial  velocity  components 
after  each  time  step,  S|u(t  +  At)  —  u{t)\  and  S|v(t  +  At)  —  v(t)|,  must  be  less  than  a 
specified  small  number.  This  approach,  employed  for  a  range  of  Mach  numbers  and 
Reynolds  numbers  with  satisfactory  results,  requires  large  computation  times. 

The  second  approach  is  to  use  the  analytically  calculated  velocity  profiles  for 
incompressible,  inviscid  flow  in  a  long,  narrow  cylindrical  tube  (Culick^)  as  start¬ 
ing  profiles  for  the  steady,  compressible,  viscous  flow  computations.  This  approach 
reduces  the  computation  time  required  to  reach  the  final  converged  steady  flow 
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configuration.  Results  given  in  Figs.  2a  and  2b  show  the  steady  normalized  axial, 
us(a:,r)/u5(i,r  =  0),  and  radial  velocity  rs(x.r)  profiles  at  different  axial  locations, 
I  =  0.025,0.5  and  1.0,  when  M  =  0.05  and  Re  —  3.10®.  respective!}'.  In  these  graphs, 
Culick®  profiles  are  indistinguishable  from  the  computed  profiles. 

The  mean  axial  flow  Mach  number  is  small,  M  =  0(10"*  -  10“^),  which  implies 
0{M^)  differences  between  the  Culick^  solution  for  ^  >>  1  and  the  computational 
result.  Nonetheless,  it  is  useful  to  compute  a  steady  state  flow  solution  for  each  Mach 
number  and  Reynolds  number.  This  prevents  introducing  unwanted  noise  into  the 
unsteady  computations. 

Although  the  flow  in  the  chamber  ,  to  a  first  order  approximation,  is  basically 
invisdd  in  character  as  a  result  of  large  Reynolds  number,  it  is  useful  to  retain  the 
radial  transport  terms  in  the  Navier-Stokes  equations.  The  steady  flow  solution  is 
obtained  faster  and,  at  the  same  time,  the  largest  important  viscous  effects  in  radial 
direction  are  responsible  for  physically  meaningful  damping,  similar  to  the  artificial 
damping  terms  that  have  been  introduced  in  some  earlier  studies^^’^^.  For  the  present 
computations  Re  number  based  on  the  characteristic  mean  axial  flow  speed  is  between 
3. 10^.10®. 

II. b.  Unsteady  Flow  Computations 

Once  a  converged  steady  flow  configuration  for  certain  M  and  Re  is  obtained,  the 
flow  is  disturbed  by  imposing  a  sinusoidally  fluctuating  component  on  the  exit  plane 
static  pressure,  a  procedure  similar  to  that  of  Vuillot  and  Avalon^®.  Therefore,  the 
the  static  boundary  condition  at  the  exit  plane  is  changed  to 

1  =  1;  p  =  1  +  Asinut  (8) 

where  u;  is  the  dimensionless  angular  frequency  and  A  is  the  amplitude  of  the  pressure 
oscillation.  The  other  boundary  conditions  are  the  same  as  those  for  steady  flow 
computations. 

The  weakly  nonlinear  theory  of  Zhao  and  Kassoy*®  and  Zhao  et  al.**  demonstrates 
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that  when  .4  =  0{M)  changes  in  the  axial  velocity  field  are  0(1),  e.g.  Au  =  0(1). 
As  a  result,  the  amplitudes  of  the  pressure  oscillations  in  the  present  computations 
are  chosen  to  be  A  =  M .  Therefore,  nonlinear  processes  affect  the  evolution  of  the 
unsteady  flow  field. 

With  reference  to  the  analytical  study  by  Zhao  et  al.*®,  it  may  be  shown  that  the 
resonant  frequencies  of  the  planar  part  of  the  acoustical  phenomenan  described  in  this 

study  are  =  (n  —  1/2)  tt  where  n  =  1,2, .  Unsteady  computations  are  carried 

for  several  different  angular  frequencies  including  those  for  near-resonant  cases. 

III.  RESULTS  and  DISCUSSION 

Following  a  procedure  described  by  Lagerstrom^®,  and  similar  to  that  employed 
by  Price  and  Flandro*®,  Zhao  and  Kassoy*^  and  Zhao  et  al.*®,  the  total  unsteady  axial 
flow  speed  may  be  divided  into  three  parts 

u(x,  r,  i  =  us{x,  r)  A  up{x,  t)  A  uv{x,  r,  t)  (9) 

where  us  denotes  the  steady  flow  field  which  is  known  as  an  initial  condition  for 
unsteady  computations,  up  represents  the  irrotational  planar  part  of  the  flow  field 
which  can  be  computed  as  the  difference  between  the  unsteady  axial  speed  and  the 
steady  axi2d  speed  on  the  centerline  of  the  tube.  The  remaining  part  uv  is  defined  as 
the  vortical  (rotational,  nonplanar)  part  of  the  unsteady  axial  flow  speed.  It  is  used  to 
describe  the  generation  and  evolution  of  the  nonlinear  unsteady  vorticity  field  in  the 
cylinder.  Following  the  aisymptotic  analysis  described  by  Zhao  et  al.*®,  one  can  show 
that  the  Vortical  part  of  the  unsteady  axial  flow  speed  uy  varnishes  at  the  centerline 
at  aiU  times. 

It  is  noted  that,  for  all  cases  to  be  discussed  below,  the  pressure  solution  has 
been  found  to  be  purely  planar  (i-dependent),  with  no  detectable  '  '•ansverse  (radial) 
variation  within  computational  accuracy.  For  example,  Fig.  3  shows  the  axial  pressure 
variation  at  three  radial  locations  (r  =  0.,0.5  and  1)  for  t  =  29.37,  M  =  0.1,  Re  = 
3.10^,  u;  =  1  and  A  =  0.1  .  In  this  figure,  numerical  pressure  values  ,  at  z  =  1/2  ,  for 
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three  radial  locations,  r  =  0,0.5  and  1  are  0.97071351,  0.97070410  and  0.97069881, 
respectively.  This  is  a  consequence  of  the  large  aspect  ratio,  S  =  20. 

Figures  4a-4c  show  the  radial  variation  of  the  unsteady  axial  vortical  flow  speed 
at  midchamber  (i  =  1/2)  at  three  time  values  after  the  planar  pressure  disturbance 
is  initiated  at  the  exit  plane.  The  flow  parameters  are  M  =  0.1  and  Re  =  3.10^.  The 
corresponding  injection  Mach  number  Mi  =  M/S  =  0.005.  The  distu-bance  frequency 
is  u;  =  1.0  ,  a  non-resonant  frequency  smaller  than  the  first  natural  frequency  of  the 
tube,  0)1  =  r/2. 

In  the  first  of  this  sequence  of  graphs,  Fig.  4a,  o  e  observes  a  strong  radial  velocity 
gradient  extending  out  about  0.15  units  from  the  wall  at  t  =  1.48.  This  is  the 
approximate  radial  distance  travelled  by  the  injected  fluid  during  the  time  interval 
t  =  1.48  to  be  discussed  at  the  end  of  this  section. 

Figure  4b  shows  ’.hat  the  unsteady  vortical  axial  velocity  field  extends  out  to  0.85 
radial  units  from  the  injecting  wall  when  t  =  11.81  .  At  time  t  =  29.37  the  rotational 
flow  field  as  seen  in  Fig.  4c  is  all  over  the  entire  chamber. 

The  spatial  distribution  of  the  vortical  part  of  the  unsteady  axial  flow  velocity  at 
each  time  may  be  explained  in  physical  terms  by  considering  an  interaction  be’  ween 
the  steady  injected  flow  field  and  the  propagating  planar  acoustic  disturbances  orig¬ 
inated  and  sustained  at  the  exit  plane  of  the  tube.  The  motion  of  a  fluid  particle 
injected  radially  into  the  tube  from  the  upper  wall  at  a  specified  location  is  affected 
by  the  harmonic  variation  with  time  of  the  local  axial  planar  pressure  gradient.  For 
instance.  Fig.  5  shows  the  time  variation  of  the  axial  pressure  gradient, dp/ dz,  at  a 
point  where  z  =  1/2  and  r  =  0.9  for  the  case  being  discussed  above,  where  the  high 
frequency  response  in  the  first  cycle,  probably  due  to  a  start-up  process,  diseappears 
quickly.  As  a  result,  a  given  fluid  particle  emanating  from  the  wall  will  be  accelerated 
alternately  in  the  positive  and  negative  axial  directions  as  it  is  convected  toward  the 
axis  of  the  cylinder  by  the  steady  radial  flow  field.  Part  of  the  fluid  particle  response 
is  associated  with  irrotational  acoustic  effects.  The  rest  is  rotational,  resulting  from 
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vorticity  generation  at  the  wall. 

Figure  6a  shows  the  time  history  of  the  planar  part  of  the  unsteady  axial  flow 
speed  at  raidchamber  when  the  flow  parameters  are  M  =  0.1,  Re  =  3.10^,  w  =  1  and 
A  =  0.1.  The  time  history  of  the  unsteady  vortical  axiad  flow  speed  at  the  injecting 
wall,  r  =  1,  when  z  =  1/2  is  shown  in  Fig.  6b.  It  is  noted  from  this  figure  that, 
before  t  =  0.5,  there  is  no  spatial  variation  of  the  unsteady  axial  vortical  speed  with 
the  radius  r  because  the  signal  originated  at  the  exit  plane  hats  not  yet  reached  the 
midchaumber.  Once  the  signal  arrives,  unsteady  vorticity  is  generated  at  the  wall 
by  the  interaction  of  the  acoustic  disturbance  amd  the  injected  gas,  as  seen  by  the 
behavior  of  uy  in  Fig.  6b.  Subsequently,  the  unsteady  vorticity  is  convected  towards 
the  axis  of  the  tube  by  the  mean  radial  flow  field,  shown  in  Fig.  2b.  There  is  also  an 
associated  pressure  field  which  has  the  characteristics  of  a  traditional  planar  acoustic 
wave  system.  This  pressure  field,  as  demonstrated  in  the  as3rmptotic  analysis  of 
Zhao  and  Kassoy^*  and  Zhao  et  al.^®,  is  not  affected  by  the  unsteady  vorticity  field 
generation  and  evolution,  so  that  traditional  acoustic  theory  pelds  transient  pressure 
estimates  that  compare  well  with  those  found  experimentally. 

Figure  7  shows  the  radial  variation  of  the  difference  between  the  unsteady  radial 
flow  speed  v  and  the  steady  radial  flow  speed  vs  in  Fig.  2b,  at  i  =  0.5  when  t  =  29.37 
for  the  same  parameter  values  above.  The  difference  is  used  because  the  as3rmptotic 
formulation  in  Zhao  and  Kassoy^^  and  Zhao  et  al.’®  demonstrates  that  the  transient 
response  to  the  boundary  disturbance  is  0{M)  smaller  than  vs  itself.  The  maximum 
absolute  .value  of  the  difference  in  Fig.  7,  approximately  0.085,  varifies  this  prediction. 

Figures  8a-8c  show  the  spatial  oscillation  of  vortical  axial  velocity  at  midchamber 
with  respect  to  the  radius  when  t  =  1.488,11.834  and  29.54  for  a  smaller  axial  Mach 
number  M  =  0.05  (corresponding  to  the  weaker  injection,  Minj  =  0.0025)  and  for  a 
higher  Re  =  3.10®.  The  forcing  frequency  u;  =  1.0  is  the  same  as  for  the  previous 
ceise.  The  amplitude  of  the  nonresonant  exit  plane  pressure  disturbance  is  5  percent. 
The  sharply  defined  region  of  large  velocity  gradient  is  seen  in  Fig.  8a  at  0.07  units 
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from  the  wall  at  f  =  1.48  .  One  notes  that  at  t  =  11.834  the  wavelength  of  the 
spatial  oscillation  of  the  vortical  axial  velocity  field  is  smaller  than  that  for  the  case 
when  M  =  0.1.  This  is  an  expected  result  because  the  mean  radial  velocity  field  for 
M  =  0.05,  which  transports  the  fluid  particles  into  the  cylinder,  is  characterized  by 
a  relatively  lower  speed  than  that  for  the  M  =  0.1  case.  Therefore,  injected  fluid 
particles  are  carried  a  shorter  distance  away  from  the  sidewall  towards  the  axis  of  the 
chamber  in  the  same  time  interval,  compared  to  that  for  the  stronger  injection  speed 
caise,  M  =  0.1.  At  t  =  29.54  one  notes  spatial  oscillations  throughout  the  cylinder. 
It  is  also  observed  that  the  wavelength  of  the  oscillatory  structure  decreases  as  the 
centerline  is  approached.  This  occurs  basically  because  the  mean  radial  flow  speed 
vanishes  as  the  centerline  is  approached. 

Solution  resolution  requires  41  grid  points  in  the  axial  direction  and  101  grid  points 
in  the  radial  direction.  Figure  8c  shows  that  near  the  injecting  wall  one  wavelength  of 
the  spatial  oscillation  of  the  vortical  axial  velocity  is  represented  by  approximately  35 
radial  grid  points.  In  contrast,  neair  the  centerline,  where  the  wavelength  is  smaller, 
approximately  10-15  grid  points  are  available  to  resolve  the  velocity  gradients. 

The  axial  spatial  variations  in  uy  do  not  have  steep  gradients,  so  that  41  equally 
spaced  grid  points  have  been  found  to  give  adequate  resolution.  For  example,  Fig.  9 
shows  the  axial  variation  of  the  unsteady  vortical  axial  velocity,  uy,  at  3  radial  lo¬ 
cations,  r  =  l.,0.9  and  0.7  when  t  =  29.37  for  M  =  0.1,  Re  =  3.10^,  u;  =  1.  and 
A  =  0.1.  In  fact,  from  Zhao  and  Kassoy”  and  Zhao  et  al.*®  one  can  show  that  the 
linear  eigenfunctions  for  the  finite  cylinder  with  a  pressure  node  exit  plane  can  be 

written  as  sin  [(2n  -  l)/4]  2irx,  n  =  1, 2, .  Thus  even  for  n=9,  where  4.25  harmonic 

waves  are  present  in  the  cylinder,  there  arc  nearly  10  grid  points  per  wave  length. 

The  third  case  studied  is  for  a  smaller  mean  axial  flow  Mach  number  M  =  0.02 
{Minj  =  0.001),  Resmolds  number  Re  =  3.10®  and  the  forcing  frequency  u;  =  1.0. 
The  results  for  the  previous  cases,  Af  =  0.1  and  M  —  0.05,  imply  that  the  number 
of  radial  grid  points  should  be  doubled  for  this  weak  injection  case.  There  are  201 
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equadly  spaced  gnd  points  in  the  radial  direction  in  order  to  represent  the  spatial 
variation  of  unsteady  vortical  axial  velocity  accurately.  Figures  lOa-lOc  show  the 
unsteady  vortical  axial  velocity  variation  with  respect  to  the  radius  at  x  =  0.5  when 
t  =  3.00,  15.00  and  30.00,  somewhat  different  from  the  previuos  cases.  It  can  be  seen 
from  these  figures  that  axial  velocity  gradients  we  larger  than  those  for  larger  Mach 
number  cases  presented  previously,  as  predicted  by  Zhao  and  Ktissoy®*  and  Zhao  et 
al.*®  This  implies  that  the  absolute  magnitude  of  the  unsteady  vorticity  generated  at 
the  wall  is  much  larger  than  that  of  the  higher  Mach  number  flows.  This  unsteady 
vorticity  field  is  convected  away  from  the  wadi  towairds  the  center  of  the  chamber  by 
a  relatively  slower  steady  radial  velocity  component.  Therefore,  at  t  =  30.00  only 
about  60  percent  of  the  chamber  is  filled  with  the  unsteady  vorticity.  Larger  time 
computations  for  this  case  have  not  been  cairhed  out.  However,  one  can  conclude  on 
the  basis  of  the  previous  results  and  the  work  of  Zhao  and  Kassoy**  and  Zhao  et  al.** 
that  the  vorticity  field  will  spread  out  towards  the  axis  as  time  increases.  One  should 
note  that  the  convection  process  slows  down  near  the  eods  because  the  radial  velocity 
component  becomes  vanishingly  small. 

Figures  11a  and  11b  present  the  spatial  vaiiation  of  the  unsteady  vortical  axial 
velocity  with  radius  at  two  different  axial  locations  x  =  1/4  and  x  =  3/4  when  t  =  30 
for  M  =  0.02.  The  corresponding  result  at  x  =  1/2  is  given  in  Fig.  10c.  It  can  be 
seen  that  the  amplitudes  of  the  oscillations  at  a  given  radial  location  decrease  as  the 
distance  from  the  source  of  disturbance  at  x  =  1  increases.  The  local  peaks  occur  at 
the  same  radial  location  for  each  axial  position. 

Figures  12a*l2c  present  the  radial  variation  of  unsteady  vorticity  at  midchamber 
length  for  three  cawes,  M  =  0.1,0.05  and  0.02,  discussed  above,  respectively.  Corre¬ 
sponding  times  are  t  =  29.37,29.54  and  30.,  respectively.  The  unsteady  vorticity  is 
computed  from  the  following  expression 

_ [9^  _  1  d{v-vs) 

dr  6^  dx 

The  non  dimensional  vorticity  is  defined  as  Q  =  fi'/  (UrIR)  where  Cl'  is  the  dimen- 
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sional  vorticity.  As  analysis  of  the  present  numerical  results  and  the  asymptotic 
analysis  of  Zhao  et  al.^*  show,  the  main  contribution  to  the  unsteady  vorticity  is 
brought  by  the  first  term  in  (10).  It  is  seen  from  Figs.  I2a-12c  that  the  magnitude  of 
the  unsteady  vorticity  increases  with  decreasing  mean  flow  Mach  number.  One  should 
notice  that  scales  are  different  in  each  of  Figs.  12a-c.  An  accompanying  plot,  Fig.  I2d, 
shows  the  variation  of  the  unsteady  vorticity  throughout  the  cylindrical  chamber  for 
the  M  —  0.02  case  at  f  =  30,  discussed  above.  The  axial  variation  of  the  vorticity  at 
each  radial  location  is  of  the  shape  of  the  eigenfunction  sin[(2n  —  l)/4|27ri  for  n  =  1, 
mentioned  previously. 

Figures  13a,  13b  and  13c  describe  the  effects  of  forcing  frequency  on  the  oscillatory 
spatial  structure  of  the  vortical  axial  speed  at  times  close  to  <  =  17.5  when  M  =  0.1 
,Re  ~  3.10^  and  A  —  0.1.  The  chosen  frequencies  are  u;  =  1.0,  u;  =  1.5  (a  near- 
resonant  case)  and  u;  =  2.5.  The  scale  in  Fig.  13b  for  a?  =  1.5  is  different  from  the 
others.  It  can  be  seen  from  these  plots  that  the  wavelength  of  the  oscillatory  spatial 
structure  decreases  when  the  forcing  frequency  increases.  Were  this  done  for  smaller 
M  the  required  number  of  radial  grid  points  would  have  to  be  significantly  larger 
to  resolve  the  growing  number  of  spatial  oscillations  at  a  given  value  of  time.  In 
addition,  larger  amounts  of  computer  time  are  required  because  the  incompressibility 
limit  is  approached  when  M  is  lowered. 

It  is  noted  from  Figs.  13a-l3c  that  the  local  spatied  gradients  of  vortical  velocity  in¬ 
crease  for  larger  forcing  frequencies.  This  implies  that  the  magnitude  of  the  unsteady 
vorticity  field  increzises  for  higher  frequencies.  One  should  notice  the  near-resonant 
forcing  frequency  effect  on  the  increasing  magnitude  of  the  unsteady  vortical  axial 
speed. 

The  radial  location  of  the  front  of  the  spatially  oscillatory  vortical  psirt  of  the 
unsteady  axial  speed  varies  with  time.  Close  to  the  sidewall  the  front  location  may 
be  computed  approximately  firom 

TF  =  l-Mt  (11) 
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where  the  mean  radial  convection  speed  is  nearly  one  in  nondimensional  terms.  Equa¬ 
tion  (11)  yields  reasonable  results  for  rp  >  0.5,  as  shown  in  Table  I.  A  more  accurate 
expression  that  defines  the  location  of  the  front  may  be  obtained  by  employing  the 
Culick*  steady  radiad  speed,  t)s(’‘)  =  — (l/r)sin  [(7r/2)  r^]  in  the  integration  of  the 
differential  form  of  time-distance  relation,  dr'  =  Vg{r')dt',  where  a  prime  denotes  a 
dimensional  quantity.  This  yields  the  radial  location  time  variation, 

rF=  -tan-*  (e-’"‘)  (12) 

Table  I  shows  the  radial  locations  estimated  from  (11)  and  (  12  )  and  from  numer¬ 
ical  computations  for  the  three  mean  axial  flow  Mach  numbers  discussed  previously. 
There  is  an  excellent  agreement  between  estimated  radial  locations  of  the  front  from 
(12)  and  from  numerical  results.  Although  (11)  gives  good  estimates  near  the  wall,  it 
yields  unreasonable  values  for  r  >  0.5.  This  is  a  result  of  the  vanishing  steady  radial 
velocity  as  r  — *  0. 

IV.  SUMMARY  and  CONCLUSIONS 

There  is  now  a  considerable  body  of  evidence  in  support  of  the  presence  of  an 
unsteady  vorticity  distribution  within  a  physically  re2isonable  model  of  a  solid  rocket 
engine  chamber.  The  experiments  of  Brown  et  al.*®’*^,  the  quasi-analytical  modeling 
of  Flandro  and  Roach**,  Zhao  and  Kassoy**,  Zhao  et  al.**  as  well  as  the  computational 
solutions  of  Vuillot  and  Avalon*®,  Flandro  and  Roach**  and  of  the  present  work  show 
unequivocally  that  unsteady  vorticity  is  generated  near  the  cylindrical  surface  and 
is  converted  into  the  chamber  by  the  injected  fluid.  The  core  of  the  chamber  is  free 
of  vorticity  only  during  the  very  early  phases  of  the  transient  process,  prior  to  the 
arrival  of  a  well  defined  unsteady  vorticity  front. (see,  for  example,  Figs.  4a,  8a  and 
10a).  The  arrival  times  found  from  the  present  computational  solution  agree  quite 
well  with  predictions  found  from  generally  valid  concepts  developed  in  the  analytical 
problem  formulation  in  Zhao  and  Kassoy**,  and  Zhao  et  al.**.  (Table  I). 

The  computations  in  the  present  work  are  probably  more  accurate  than  those  of 
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Vuillot  and  Avalon^®  because  the  analytically  derived  scaling  is  used  to  determine 
grid  size  and  distribution.  Unlike  the  former  work,  where  grid  points  were  concen¬ 
trated  near  the  sidewall  to  resolve  an  expected  traditional  acoustic  boundary  layer, 
the  present  choices  aissure  that  the  short  wavelength  processes  in  the  radial  spatial 
distributions,  predicted  by  the  analysis,  are  adequately  resolved.  The  time-averaged 
axial  velocity  variation  with  radius  predicted  by  Vuillot  and  Avalon^®  (see  their  Fig.  7) 
is  qualitatively  similar  to  the  instantaneous  distributions  for  the  vortical  axial  velocity 
described  here. 

One  may  reasonably  conclude  from  the  results  of  the  aforementioned  studies  that 
the  classical  acoustic  stability  predictions  used  for  scientific  and  engineering  purposes 
should  be  reassessed  to  determine  the  prediction  reliability.  It  is  conceivable  that  the 
acoustic  pressure  predictions  are  reasonable  because  the  presence  of  vorticity  does  not 
impact  the  primary  acoustic  pressure  transient  found  in  the  Zhao  and  Kassoy^*  and 
Zhao  et  al.*®  theory.  However,  it  is  unlikely  that  velocity  and  shear  stress  predictions 
are  accurate,  given  the  irrotational  basis  of  the  acoustic  stability  analysis. 

It  is  also  conceivable  that  there  are  parameter  ranges  of  Mach  number  (injection 
rate)  and  Reynolds  number  for  which  vorticity  is  really  confined  to  thin,  classical, 
but  injected  acoustic  boundary  layers  like  those  discussed  by  Flandro^^,  Baum  and 
Levine^®  and  more  recently  by  Zhao  et  al.*®.  These  structures  can  appear  for  relatively 
small  injection  rates  and  low  Reynolds  numbers,  so  that  viscous  damping  of  the 
vorticity  amplitude  is  profound.  Then,  the  cylinder  core  will  contain  the  relatively 
weak  vorticity  of  the  steady  Culick^  solution  and  irrotational  acoustic  waves  driven 
by  the  boundary  forcing. 

The  amplitude  of  the  transient  vorticity  distributions  described  by  Zhao  and 
Kassoy*^  and  Zhao  et  al.*®,  and  in  the  present  work  (see  Figs.  12a-c)  are  0{M~^) 
larger  than  that  of  the  Culick^  steady  solution.  This  implies  that  there  will  be  a 
relatively  large  axial  shear  stress  on  the  cylinder  surface,  particularly  for  smaller  M 
values.  One  can  speculate  that  these  transient  shear  stresses  will  impact  the  burning 
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rate  of  an  actual  propellant  which  is  the  source  of  the  “injected”  fluid  used  in  the 
present  model.  Perhaps  there  is  a  direct  relationship  between  the  effect  of  surface 
shear  stress  transients,  predicted  in  the  present  work,  and  erosive  burning  concepts 
used  in  the  sobd  rocket  engineering  literature**. 
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TABLE  CAPTION 


Table  I:  The  radial  locations  of  the  unsteady  vortical  axial  velocity  front  at 
different  time  levels  for  M  =  0.02,0.05  and  0.1  .  The  second  and  third  columns 
present  the  estimates  from  (11)  and  (12),  respectively.  Results  in  the  leist  column 
have  been  found  from  Figs.  4,  8  and  10. 
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FIGURE  CAPTIONS 


Fig.  1:  Computational  domain  and  boundary  conditions. 


Fig.  2a:  Normalized  steady  ajcial  velocity  profiles  at  z  =  0.025,  0.5  and  I  for 
M  =  0.05  and  Re  =  3.10®  . 


Fig.  2b:  Steady  radial  velocity  profiles  at  axial  locations  z  =  0.025  (solid  line),  0.5 
(dotted  line)  and  1  (daished  line)  for  M  —  0.05  and  Re  =  3.10®. 


Fig.  3:  The  axial  unsteady  pressure  variation  at  r  =  0,  0.5  and  1  when  t  =  29.37 
for  Af  =  0.1,  Re  =  3.10^,  u;  =  1  and  A  =  M. 


Fig.  4a:  The  radial  variation  of  the  unsteady  axial  flow  speed, av,  at  z  =  0.5  when 
t  =  1.48  for  M  =  0.1,  Re  =  3.10^,  u;  =  1  and  A  =  M. 


Fig.  4b:  As  Fig.  4a  but  when  t  =  11.81. 


Fig.  4c:  As  Fig.  4a  but  when  t  =  29.37. 


Fig.  5:  The  time  history  of  axial  pressure  gradient,  dpfdx,  at  z  =  0.5,  r  =  0.9  for 
M  =  0.1,  Re  =  3.10^  u;  =  1  and  A  =  M. 
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Fig.  6a:  The  time  history  of  the  planar  part  of  the  unsteady  axial  flow  speed,  up 
at  z  =  0.5  for  A/  =  0.1,  Re  =  3.10^,  u;  =  1  and  A  =  M. 


Fig.  6b:  The  time  history  of  the  unsteady  vortical  axizd  speed,  uv  at  z  =  0.5, 
r  =  1  for  the  same  flow  parameters  as  those  in  Fig.  6a. 


Fig.  7:  The  radial  variation  of  (u  — 115)  at  z  =  0.5  when  t  =  29.37  for  the  same  flow 
parameters  as  those  in  Fig.  4c. 


Fig,  8a:  The  radial  variation  of  uv  at  z  =  0.5  when  t  =  1.488  for  M  —  0.05, 
Re  -  3.10®,  O'  =  1  and  A  =  M. 


Fig.  8b:  As  Fig.  8a  but  for  t  =  11.834. 


Fig.  8c:  As  Fig.  8a  but  for  t  —  29.54. 


Fig.  9:  The  eixial  variation  of  uv  at  r  =  1  (solid  line),  r  =  0.9  (dotted  line)  and 
r  =  0.7  (  dashed  line)  when  t  =  29.37  for  M  =  0.1,  Re  =  3.10^,  u;  =  1  and  A  =  M. 


Fig.  10a:  The  radial  variation  of  uv  at  i  =  0.5  when  f  =  3  for  Af  =  0.02, 
Re  —  3.10®,  u;  =  1  and  A  =  M. 


Fig.  10b:  As  Fig.  10a  but  when  t  =  15. 
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Fig.  10c:  As  Fig  10a  but  when  t  =  30. 

Fig.  11a:  As  Fig.  10c  but  at  i  =  0.25. 

Fig.  lib:  As  Fig.  11a  but  at  x  =  0.75. 

Fig.  12a:  The  radial  variation  of  unsteady  vorticity,  fi,  defined  in  (10)  at  i  =  0.5 
when  t  =  29.37  for  A/  =  0.1  . 

Fig.  12b:  As  Fig.  12a  but  for  M  —  0  05,  when  t  =  29.54. 

Fig.  12c:  As  Fig.  12a  but  for  M  =  0.02,  when  t  =  30. 

Fig.  I2d:  Spatial  unsteady  ’orticity  variation  Ct  throughout  the  cylinder  at  t  =  30 
for  M  =  0.02. 

Fig.  13a:  The  radial  variation  of  uv  x  =  0.5  when  t  =  17.71  for  M  =  0.1, 

Re  =  3.10*;  A  =  M  and  u;  =  1. 

Fig.  13b:  As  Fig.  13a  but  at  t  =  17.38  and  for  uf  =  1.5. 

Fig.  13c:  As  Fig.  13a  but  at  t  =  17.77  and  for  u/  =  2.5. 
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